Context Navigation

physiq.F @ 1015

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

Ajoute de initialisation a zero pour 2 variables qui ne sont pas calcule dans le cas de Tietke.
Copure des lignes trop longues
Correction de nom pour 2 variables

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

File size: 111.0 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 zustrhi(klon), zvstrhi(klon)
946	REAL aam, torsfc
947	cIM 141004 END
948	cIM 190504 BEG
949	INTEGER ij, imp1jmp1
950	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
951	cym A voir plus tard
952	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
953	REAL padyn(iim+1,jjmp1,klev+1)
954	REAL dudyn(iim+1,jjmp1,klev)
955	REAL rlatdyn(iim+1,jjmp1)
956	cIM 190504 END
957	LOGICAL ok_msk
958	REAL msk(klon)
959	cIM
960	REAL airetot, pi
961	cym A voir plus tard
962	cym REAL zm_wo(jjmp1, klev)
963	cIM AMIP2 END
964	c
965	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
966	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
967	c#ifdef histmthNMC
968	cym A voir plus tard !!!!
969	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
970	REAL zx_tmp_fiNC(klon,nlevSTD)
971	c#endif
972	REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
973	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
974	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
975	c
976	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
977	INTEGER nid_ctesGCM
978	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
979	SAVE nid_ctesGCM
980	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
981	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
982	c
983	cIM 280405 BEG
984	INTEGER nid_bilKPins, nid_bilKPave
985	SAVE nid_bilKPins, nid_bilKPave
986	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
987	c
988	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
989	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
990	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
991	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
992	c
993	cIM 280405 END
994	c
995	INTEGER nhori, nvert, nvert1, nvert3
996	REAL zsto, zsto1, zsto2
997	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
998	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
999	REAL zout_isccp(napisccp)
1000	SAVE zcals, zcalh, zoutj, zout_isccp
1001	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1002
1003	real zjulian
1004	save zjulian
1005	c$OMP THREADPRIVATE(zjulian)
1006
1007	character*20 modname
1008	character*80 abort_message
1009	logical ok_sync
1010	real date0
1011	integer idayref
1012
1013	C essai writephys
1014	integer fid_day, fid_mth, fid_ins
1015	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1016	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1017	parameter (prof2d_on = 1, prof3d_on = 2,
1018	. prof2d_av = 3, prof3d_av = 4)
1019	character*30 nom_fichier
1020	character*10 varname
1021	character*40 vartitle
1022	character*20 varunits
1023	C Variables liees au bilan d'energie et d'enthalpi
1024	REAL ztsol(klon)
1025	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1026	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1027	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1028	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1029	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1030	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1031	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1032	REAL d_h_vcol_phy
1033	REAL fs_bound, fq_bound
1034	SAVE d_h_vcol_phy
1035	c$OMP THREADPRIVATE(d_h_vcol_phy)
1036	REAL zero_v(klon)
1037	CHARACTER*15 ztit
1038	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1039	SAVE ip_ebil
1040	DATA ip_ebil/0/
1041	c$OMP THREADPRIVATE(ip_ebil)
1042	INTEGER if_ebil ! level for energy conserv. dignostics
1043	SAVE if_ebil
1044	c$OMP THREADPRIVATE(if_ebil)
1045	c+jld ec_conser
1046	REAL ZRCPD
1047	c-jld ec_conser
1048	REAL t2m(klon,nbsrf) ! temperature a 2m
1049	REAL q2m(klon,nbsrf) ! humidite a 2m
1050
1051	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1052	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1053	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1054	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1055	CHARACTER*40 tinst, tave, typeval
1056	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
1057	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
1058
1059	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1060
1061	REAL re(klon, klev) ! Cloud droplet effective radius
1062	REAL fl(klon, klev) ! denominator of re
1063
1064	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1065
1066	! Aerosol optical properties
1067
1068	! Aerosol optical properties by INCA model
1069	CHARACTER*4 :: rfname(9)
1070	REAL aerindex(klon) ! POLDER aerosol index
1071
1072	! Parameters
1073	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1074	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1075	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1076	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1077	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1078	! false : lecture des aerosol dans un fichier
1079	c$OMP THREADPRIVATE(aerosol_couple)
1080
1081	c
1082	c Declaration des constantes et des fonctions thermodynamiques
1083	c
1084	LOGICAL,SAVE :: first=.true.
1085	c$OMP THREADPRIVATE(first)
1086	#include "YOMCST.h"
1087	#include "YOETHF.h"
1088	#include "FCTTRE.h"
1089	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1090	#include "conema3.h"
1091	#include "fisrtilp.h"
1092	#include "nuage.h"
1093	#include "compbl.h"
1094	cIM 100106 END : pouvoir sortir les ctes de la physique
1095	c
1096	c======================================================================
1097	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1098	! Utilise notamment en 1D mais peut etre active egalement en 3D
1099	! en imposant la valeur de igout.
1100	c======================================================================
1101
1102	if (prt_level.ge.1) then
1103	igout=klon/2+1/klon
1104	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1105	write(lunout,*)
1106	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
1107	write(lunout,*)
1108	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
1109
1110	write(lunout,*) 'papers, play, phi, u, v, t, omega'
1111	do k=1,nlev
1112	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1113	s u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
1114	enddo
1115	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1116	do k=1,nlev
1117	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1118	enddo
1119	endif
1120
1121	c======================================================================
1122
1123	cym => necessaire pour iflag_con != 2
1124	pmfd(:,:) = 0.
1125	pen_u(:,:) = 0.
1126	pen_d(:,:) = 0.
1127	pde_d(:,:) = 0.
1128	pde_u(:,:) = 0.
1129	aam=0.
1130	torsfc=0.
1131
1132	if (first) then
1133
1134	cCR:nvelles variables convection/poches froides
1135
1136	print*, '================================================='
1137	print*, 'Allocation des variables locales et sauvegardees'
1138	call phys_local_var_init
1139	call phys_state_var_init
1140	print*, '================================================='
1141
1142	cIM beg
1143	dnwd0=0.0
1144	ftd=0.0
1145	fqd=0.0
1146	cin=0.
1147	cym Attention pbase pas initialise dans concvl !!!!
1148	pbase=0
1149	paire_ter(:)=0.
1150	cIM 180608
1151	c pmflxr=0.
1152	c pmflxs=0.
1153	first=.false.
1154
1155	endif ! fisrt
1156
1157	modname = 'physiq'
1158	cIM
1159	IF (ip_ebil_phy.ge.1) THEN
1160	DO i=1,klon
1161	zero_v(i)=0.
1162	END DO
1163	END IF
1164	ok_sync=.TRUE.
1165	IF (nqmax .LT. 2) THEN
1166	abort_message = 'eaux vapeur et liquide sont indispensables'
1167	CALL abort_gcm (modname,abort_message,1)
1168	ENDIF
1169	IF (debut) THEN
1170	CALL suphel ! initialiser constantes et parametres phys.
1171	ENDIF
1172
1173	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1174
1175
1176	c======================================================================
1177	xjour = rjourvrai
1178	c
1179	c Si c'est le debut, il faut initialiser plusieurs choses
1180	c ********
1181	c
1182	IF (debut) THEN
1183	C
1184	!rv
1185	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1186	cde la convection a partir des caracteristiques du thermique
1187	wght_th(:,:)=1.
1188	lalim_conv(:)=1
1189	cRC
1190	u10m(:,:)=0.
1191	v10m(:,:)=0.
1192	piz_ae(:,:,:)=0.
1193	tau_ae(:,:,:)=0.
1194	cg_ae(:,:,:)=0.
1195	rain_con(:)=0.
1196	snow_con(:)=0.
1197	bl95_b0=0.
1198	bl95_b1=0.
1199	topswai(:)=0.
1200	topswad(:)=0.
1201	solswai(:)=0.
1202	solswad(:)=0.
1203
1204	IF (config_inca /= 'none') THEN
1205	tau_inca(:,:,:,:) = 0.
1206	piz_inca(:,:,:,:) = 0.
1207	cg_inca(:,:,:,:) = 0.
1208	ccm(:,:,:) = 0.
1209	topswai_inca(:) = 0.
1210	topswad_inca(:) = 0.
1211	topswad0_inca(:) = 0.
1212	topsw_inca(:,:) = 0.
1213	topsw0_inca(:,:) = 0.
1214	solswai_inca(:) = 0.
1215	solswad_inca(:) = 0.
1216	solswad0_inca(:) = 0.
1217	solsw_inca(:,:) = 0.
1218	solsw0_inca(:,:) = 0.
1219	END IF
1220
1221	rnebcon0(:,:) = 0.0
1222	clwcon0(:,:) = 0.0
1223	rnebcon(:,:) = 0.0
1224	clwcon(:,:) = 0.0
1225
1226	cIM
1227	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1228	c
1229	c appel a la lecture du run.def physique
1230	c
1231	call conf_phys(ok_journe, ok_mensuel,
1232	. ok_instan, ok_hf,
1233	. solarlong0,seuil_inversion,
1234	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1235	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
1236	. ok_ade, ok_aie, aerosol_couple,
1237	. bl95_b0, bl95_b1,
1238	. iflag_thermals,nsplit_thermals,tau_thermals,
1239	cnv flags pour la convection et les poches froides
1240	. iflag_coupl,iflag_clos,iflag_wake)
1241
1242	print*,'iflag_coupl,iflag_clos,iflag_wake',
1243	. iflag_coupl,iflag_clos,iflag_wake
1244	print*,'CYCLE_DIURNE', cycle_diurne
1245
1246	c
1247	c
1248	c Initialiser les compteurs:
1249	c
1250	itap = 0
1251	itaprad = 0
1252
1253	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1254	!! Un petit travail à faire ici.
1255	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1256
1257	if (iflag_pbl>1) then
1258	PRINT*, "Using method MELLOR&YAMADA"
1259	endif
1260
1261	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1262	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1263	! dyn3d
1264	! Attention : la version precedente n'etait pas tres propre.
1265	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1266	! pour obtenir le meme resultat.
1267	dtime=pdtphys
1268	radpas = NINT( 86400./dtime/nbapp_rad)
1269	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1270
1271	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1272	cIM begin
1273	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1274	$,ratqs(1,1)
1275	cIM end
1276
1277
1278
1279	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1280	c
1281	C on remet le calendrier a zero
1282	c
1283	IF (raz_date .eq. 1) THEN
1284	itau_phy = 0
1285	ENDIF
1286
1287	cIM cf. AM 081204 BEG
1288	PRINT*,'cycle_diurne3 =',cycle_diurne
1289	cIM cf. AM 081204 END
1290	c
1291	CALL printflag( tabcntr0,radpas,ok_journe,
1292	, ok_instan, ok_region )
1293	c
1294	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1295	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1296	. pdtphys
1297	abort_message='Pas physique n est pas correct '
1298	! call abort_gcm(modname,abort_message,1)
1299	dtime=pdtphys
1300	ENDIF
1301	IF (nlon .NE. klon) THEN
1302	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1303	. klon
1304	abort_message='nlon et klon ne sont pas coherents'
1305	call abort_gcm(modname,abort_message,1)
1306	ENDIF
1307	IF (nlev .NE. klev) THEN
1308	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1309	. klev
1310	abort_message='nlev et klev ne sont pas coherents'
1311	call abort_gcm(modname,abort_message,1)
1312	ENDIF
1313	c
1314	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1315	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1316	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1317	abort_message='Nbre d appels au rayonnement insuffisant'
1318	call abort_gcm(modname,abort_message,1)
1319	ENDIF
1320	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1321	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1322	. ok_cvl
1323	c
1324	cKE43
1325	c Initialisation pour la convection de K.E. (sb):
1326	IF (iflag_con.GE.3) THEN
1327
1328	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1329	WRITE(lunout,*)
1330	. "On va utiliser le melange convectif des traceurs qui"
1331	WRITE(lunout,*)"est calcule dans convect4.3"
1332	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1333
1334	DO i = 1, klon
1335	ema_cbmf(i) = 0.
1336	ema_pcb(i) = 0.
1337	ema_pct(i) = 0.
1338	ema_workcbmf(i) = 0.
1339	ENDDO
1340	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1341	DO i = 1, klon
1342	ibas_con(i) = 1
1343	itop_con(i) = 1
1344	ENDDO
1345	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1346	c===============================================================================
1347	cCR:04.12.07: initialisations poches froides
1348	c Controle de ALE et ALP pour la fermeture convective (jyg)
1349	if (iflag_wake.eq.1) then
1350	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1351	s ,alp_bl_prescr, ale_bl_prescr)
1352	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1353	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1354	endif
1355
1356	do i = 1,klon
1357	Ale_bl(i)=0.
1358	Alp_bl(i)=0.
1359	enddo
1360
1361	c================================================================================
1362
1363	ENDIF
1364
1365	rugoro=0.
1366	c34EK
1367	IF (ok_orodr) THEN
1368
1369	rugoro=0.
1370
1371	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1372	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1373	! justement quand ok_orodr = false.
1374	! ce rugoro est utilise par la couche limite et fait double emploi
1375	! avec les paramétrisations spécifiques de Francois Lott.
1376	! DO i=1,klon
1377	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1378	! ENDDO
1379	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1380	IF (ok_strato) THEN
1381	CALL SUGWD_strato(klon,klev,paprs,pplay)
1382	ELSE
1383	CALL SUGWD(klon,klev,paprs,pplay)
1384	ENDIF
1385
1386	DO i=1,klon
1387	zuthe(i)=0.
1388	zvthe(i)=0.
1389	if(zstd(i).gt.10.)then
1390	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1391	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1392	endif
1393	ENDDO
1394	ENDIF
1395	c
1396	c
1397	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1398	WRITE(lunout,*)'La frequence de lecture surface est de ',
1399	. lmt_pas
1400	c
1401	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1402	c IF (ok_mensuel) THEN
1403	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1404	c . ecrit_mth
1405	c ENDIF
1406	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1407	c IF (ok_journe) THEN
1408	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1409	c . ecrit_day
1410	c ENDIF
1411	cIM 130904 BEG
1412	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1413	cIM 170305
1414	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1415	cIM 230305
1416	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1417	c
1418	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1419	c
1420	cIM200505 IF (ok_journe) THEN
1421	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1422	cIM200505 . ecrit_hf
1423	cIM200505 ENDIF
1424	cIM 130904 END
1425	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1426	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1427	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1428	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1429	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1430	cIM200505 IF (ok_instan) THEN
1431	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1432	cIM200505 . ecrit_ins
1433	cIM200505 ENDIF
1434	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1435	cIM200505 IF (ok_region) THEN
1436	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1437	cIM200505 . ecrit_reg
1438	cIM200505 ENDIF
1439	cIM 030306 BEG
1440	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1441	cIM : ne pas modifier ecrit_hf2mth
1442	c
1443	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1444	ecrit_hf2mth = ecrit_mth/ecrit_hf
1445	c ecrit_ins en secondes, chaque pas de temps de la physique
1446	ecrit_ins = dtime
1447	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1448	ecrit_hf = ecrit_hf * un_jour
1449	!IM
1450	IF(ecrit_day.LE.1.) THEN
1451	ecrit_day = ecrit_day * un_jour !en secondes
1452	ENDIF
1453	!IM
1454	ecrit_mth = ecrit_mth * un_jour
1455	ecrit_reg = ecrit_reg * un_jour
1456	ecrit_tra = ecrit_tra * un_jour
1457	ecrit_ISCCP = ecrit_ISCCP * un_jour
1458	c
1459	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1460	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1461	. ecrit_hf2mth
1462	cIM 030306 END
1463
1464	capemaxcels = 't_max(X)'
1465	t2mincels = 't_min(X)'
1466	t2maxcels = 't_max(X)'
1467	tinst = 'inst(X)'
1468	tave = 'ave(X)'
1469	cIM cf. AM 081204 BEG
1470	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1471	cIM cf. AM 081204 END
1472	c
1473	c=============================================================
1474	c Initialisation des sorties
1475	c=============================================================
1476
1477	#ifdef CPP_IOIPSL
1478
1479	c Commente par abderrahmane 11 2 08
1480	c#ifdef histhf
1481	c#include "ini_histhf.h"
1482	c#endif
1483
1484	c#ifdef histday
1485	c#include "ini_histday.h"
1486	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1487	c#include "ini_bilKP_ins.h"
1488	c#include "ini_bilKP_ave.h"
1489	c#endif
1490
1491	c#ifdef histmth
1492	c#include "ini_histmth.h"
1493	c#endif
1494
1495	c#ifdef histins
1496	c#include "ini_histins.h"
1497	c#endif
1498
1499	c$OMP MASTER
1500	call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta,
1501	& ctetaSTD,dtime,presnivs,ok_veget,
1502	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1503	& ok_hf,ok_instan,nid_files)
1504	c$OMP END MASTER
1505	c$OMP BARRIER
1506
1507	#ifdef histISCCP
1508	#include "ini_histISCCP.h"
1509	#endif
1510
1511	#ifdef histmthNMC
1512	#include "ini_histmthNMC.h"
1513	#endif
1514
1515	#include "ini_histday_seri.h"
1516
1517	#include "ini_paramLMDZ_phy.h"
1518
1519	#endif
1520
1521	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1522	date0 = zjulian
1523	C date0 = day_ini
1524	WRITE(,) 'physiq date0 : ',date0
1525	c
1526	c
1527	c
1528	c Prescrire l'ozone dans l'atmosphere
1529	c
1530	c
1531	cc DO i = 1, klon
1532	cc DO k = 1, klev
1533	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1534	cc ENDDO
1535	cc ENDDO
1536	c
1537	IF (config_inca /= 'none') THEN
1538	#ifdef INCA
1539	CALL VTe(VTphysiq)
1540	CALL VTb(VTinca)
1541	iii = MOD(NINT(xjour),360)
1542	calday = FLOAT(iii) + gmtime
1543	WRITE(lunout,*) 'initial time ', xjour, calday
1544
1545	CALL chemini(
1546	$ rg,
1547	$ ra,
1548	$ airephy,
1549	$ rlat,
1550	$ rlon,
1551	$ presnivs,
1552	$ calday,
1553	$ klon,
1554	$ nqmax,
1555	$ pdtphys,
1556	$ annee_ref,
1557	$ day_ini)
1558
1559	CALL VTe(VTinca)
1560	CALL VTb(VTphysiq)
1561	#endif
1562	END IF
1563	c
1564	c
1565	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1566	! Nouvelle initialisation pour le rayonnement RRTM
1567	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1568
1569	call iniradia(klon,klev,paprs(1,1:klev+1))
1570
1571	ENDIF
1572	!
1573	! ************** Fin de IF ( debut ) *************
1574	!
1575	!
1576	! Incrementer le compteur de la physique
1577	!
1578	itap = itap + 1
1579	julien = MOD(NINT(xjour),360)
1580	if (julien .eq. 0) julien = 360
1581
1582	!
1583	! Update fraction of the sub-surfaces (pctsrf) and
1584	! initialize, where a new fraction has appeared, all variables depending
1585	! on the surface fraction.
1586	!
1587	CALL change_srf_frac(itap, dtime, julien,
1588	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1589
1590
1591	! Tendances bidons pour les processus qui n'affectent pas certaines
1592	! variables.
1593	du0(:,:)=0.
1594	dv0(:,:)=0.
1595	dq0(:,:)=0.
1596	dql0(:,:)=0.
1597	c
1598	c Mettre a zero des variables de sortie (pour securite)
1599	c
1600	DO i = 1, klon
1601	d_ps(i) = 0.0
1602	ENDDO
1603	DO k = 1, klev
1604	DO i = 1, klon
1605	d_t(i,k) = 0.0
1606	d_u(i,k) = 0.0
1607	d_v(i,k) = 0.0
1608	ENDDO
1609	ENDDO
1610	DO iq = 1, nqmax
1611	DO k = 1, klev
1612	DO i = 1, klon
1613	d_qx(i,k,iq) = 0.0
1614	ENDDO
1615	ENDDO
1616	ENDDO
1617	da(:,:)=0.
1618	mp(:,:)=0.
1619	phi(:,:,:)=0.
1620	c
1621	c Ne pas affecter les valeurs entrees de u, v, h, et q
1622	c
1623	DO k = 1, klev
1624	DO i = 1, klon
1625	t_seri(i,k) = t(i,k)
1626	u_seri(i,k) = u(i,k)
1627	v_seri(i,k) = v(i,k)
1628	q_seri(i,k) = qx(i,k,ivap)
1629	ql_seri(i,k) = qx(i,k,iliq)
1630	qs_seri(i,k) = 0.
1631	ENDDO
1632	ENDDO
1633	IF (nqmax.GE.3) THEN
1634	DO iq = 3, nqmax
1635	DO k = 1, klev
1636	DO i = 1, klon
1637	tr_seri(i,k,iq-2) = qx(i,k,iq)
1638	ENDDO
1639	ENDDO
1640	ENDDO
1641	ELSE
1642	DO k = 1, klev
1643	DO i = 1, klon
1644	tr_seri(i,k,1) = 0.0
1645	ENDDO
1646	ENDDO
1647	ENDIF
1648	C
1649	DO i = 1, klon
1650	ztsol(i) = 0.
1651	ENDDO
1652	DO nsrf = 1, nbsrf
1653	DO i = 1, klon
1654	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1655	ENDDO
1656	ENDDO
1657	cIM
1658	IF (ip_ebil_phy.ge.1) THEN
1659	ztit='after dynamic'
1660	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1661	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1662	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1663	C Comme les tendances de la physique sont ajoute dans la dynamique,
1664	C on devrait avoir que la variation d'entalpie par la dynamique
1665	C est egale a la variation de la physique au pas de temps precedent.
1666	C Donc la somme de ces 2 variations devrait etre nulle.
1667	call diagphy(airephy,ztit,ip_ebil_phy
1668	e , zero_v, zero_v, zero_v, zero_v, zero_v
1669	e , zero_v, zero_v, zero_v, ztsol
1670	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1671	s , fs_bound, fq_bound )
1672	END IF
1673
1674	c Diagnostiquer la tendance dynamique
1675	c
1676	IF (ancien_ok) THEN
1677	DO k = 1, klev
1678	DO i = 1, klon
1679	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1680	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1681	ENDDO
1682	ENDDO
1683	ELSE
1684	DO k = 1, klev
1685	DO i = 1, klon
1686	d_t_dyn(i,k) = 0.0
1687	d_q_dyn(i,k) = 0.0
1688	ENDDO
1689	ENDDO
1690	ancien_ok = .TRUE.
1691	ENDIF
1692	c
1693	c Ajouter le geopotentiel du sol:
1694	c
1695	DO k = 1, klev
1696	DO i = 1, klon
1697	zphi(i,k) = pphi(i,k) + pphis(i)
1698	ENDDO
1699	ENDDO
1700	c
1701	c Verifier les temperatures
1702	c
1703	cIM BEG
1704	IF (check) THEN
1705	amn=MIN(ftsol(1,is_ter),1000.)
1706	amx=MAX(ftsol(1,is_ter),-1000.)
1707	DO i=2, klon
1708	amn=MIN(ftsol(i,is_ter),amn)
1709	amx=MAX(ftsol(i,is_ter),amx)
1710	ENDDO
1711	c
1712	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1713	ENDIF !(check) THEN
1714	cIM END
1715	c
1716	CALL hgardfou(t_seri,ftsol,'debutphy')
1717	c
1718	cIM BEG
1719	IF (check) THEN
1720	amn=MIN(ftsol(1,is_ter),1000.)
1721	amx=MAX(ftsol(1,is_ter),-1000.)
1722	DO i=2, klon
1723	amn=MIN(ftsol(i,is_ter),amn)
1724	amx=MAX(ftsol(i,is_ter),amx)
1725	ENDDO
1726	c
1727	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1728	ENDIF !(check) THEN
1729	cIM END
1730	c
1731	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1732	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1733	c
1734	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1735	WRITE(lunout,*)' PHYS cond julien ',julien
1736	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1737	ENDIF
1738	c
1739	c Re-evaporer l'eau liquide nuageuse
1740	c
1741	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1742	DO i = 1, klon
1743	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1744	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1745	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1746	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1747	zb = MAX(0.0,ql_seri(i,k))
1748	za = - MAX(0.0,ql_seri(i,k))
1749	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1750	t_seri(i,k) = t_seri(i,k) + za
1751	q_seri(i,k) = q_seri(i,k) + zb
1752	ql_seri(i,k) = 0.0
1753	d_t_eva(i,k) = za
1754	d_q_eva(i,k) = zb
1755	ENDDO
1756	ENDDO
1757	cIM
1758	IF (ip_ebil_phy.ge.2) THEN
1759	ztit='after reevap'
1760	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1761	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1762	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1763	call diagphy(airephy,ztit,ip_ebil_phy
1764	e , zero_v, zero_v, zero_v, zero_v, zero_v
1765	e , zero_v, zero_v, zero_v, ztsol
1766	e , d_h_vcol, d_qt, d_ec
1767	s , fs_bound, fq_bound )
1768	C
1769	END IF
1770
1771	c
1772	c=========================================================================
1773	! Calculs de l'orbite.
1774	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1775	! doit donc etre placé avant radlwsw et pbl_surface
1776
1777	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1778	! solarlong0
1779
1780	if (solarlong0<-999.) then
1781	CALL orbite(FLOAT(julien),zlongi,dist)
1782	else
1783	zlongi=solarlong0 ! longitude solaire vraie
1784	dist=1. ! distance au soleil / moyenne
1785	endif
1786
1787	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
1788
1789	! Avec ou sans cycle diurne
1790	IF (cycle_diurne) THEN
1791	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1792	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1793	ELSE
1794	rmu0 = -999.999
1795	ENDIF
1796
1797	if (mydebug) then
1798	call writefield_phy('u_seri',u_seri,llm)
1799	call writefield_phy('v_seri',v_seri,llm)
1800	call writefield_phy('t_seri',t_seri,llm)
1801	call writefield_phy('q_seri',q_seri,llm)
1802	endif
1803
1804	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1805	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1806	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1807	c turbulent du couche limit.
1808	c
1809	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1810	c ecriture des fihiers hist_XXXX.nc, ces sont :
1811	c qsol, zq2m, s_pblh, s_lcl,
1812	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1813	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1814	c zxrugs, zu10m, zv10m, fder,
1815	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1816	c frugs, agesno, fsollw, fsolsw,
1817	c d_ts, fevap, fluxlat, t2m,
1818	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1819	c
1820	c Certains ne sont pas utiliser du tout :
1821	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1822	c
1823
1824	CALL pbl_surface(
1825	e dtime, date0, itap, julien,
1826	e debut, lafin,
1827	e rlon, rlat, rugoro, rmu0,
1828	e rain_fall, snow_fall, solsw, sollw,
1829	e t_seri, q_seri, u_seri, v_seri,
1830	e pplay, paprs, pctsrf,
1831	+ ftsol, falb1, falb2, u10m, v10m,
1832	s sollwdown, cdragh, cdragm, yu1, yv1,
1833	s albsol1, albsol2, sens, evap,
1834	s zxtsol, zxfluxlat, zt2m, qsat2m,
1835	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1836	s ycoefh, slab_wfbils,
1837	d qsol, zq2m, s_pblh, s_lcl,
1838	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1839	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1840	d zxrugs, zu10m, zv10m, fder,
1841	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1842	d frugs, agesno, fsollw, fsolsw,
1843	d d_ts, fevap, fluxlat, t2m,
1844	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1845	- dsens, devap, zxsnow,
1846	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1847
1848
1849	!-----------------------------------------------------------------------------------------
1850	! ajout des tendances de la diffusion turbulente
1851	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1852	!-----------------------------------------------------------------------------------------
1853
1854	if (mydebug) then
1855	call writefield_phy('u_seri',u_seri,llm)
1856	call writefield_phy('v_seri',v_seri,llm)
1857	call writefield_phy('t_seri',t_seri,llm)
1858	call writefield_phy('q_seri',q_seri,llm)
1859	endif
1860
1861
1862	IF (ip_ebil_phy.ge.2) THEN
1863	ztit='after surface_main'
1864	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1865	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1866	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1867	call diagphy(airephy,ztit,ip_ebil_phy
1868	e , zero_v, zero_v, zero_v, zero_v, sens
1869	e , evap , zero_v, zero_v, ztsol
1870	e , d_h_vcol, d_qt, d_ec
1871	s , fs_bound, fq_bound )
1872	END IF
1873
1874	c =================================================================== c
1875	c Calcul de Qsat
1876
1877	DO k = 1, klev
1878	DO i = 1, klon
1879	zx_t = t_seri(i,k)
1880	IF (thermcep) THEN
1881	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1882	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1883	zx_qs = MIN(0.5,zx_qs)
1884	zcor = 1./(1.-retv*zx_qs)
1885	zx_qs = zx_qs*zcor
1886	ELSE
1887	IF (zx_t.LT.t_coup) THEN
1888	zx_qs = qsats(zx_t)/pplay(i,k)
1889	ELSE
1890	zx_qs = qsatl(zx_t)/pplay(i,k)
1891	ENDIF
1892	ENDIF
1893	zqsat(i,k)=zx_qs
1894	ENDDO
1895	ENDDO
1896
1897	if (prt_level.ge.1) then
1898	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1899	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1900	endif
1901	c
1902	c Appeler la convection (au choix)
1903	c
1904	DO k = 1, klev
1905	DO i = 1, klon
1906	conv_q(i,k) = d_q_dyn(i,k)
1907	. + d_q_vdf(i,k)/dtime
1908	conv_t(i,k) = d_t_dyn(i,k)
1909	. + d_t_vdf(i,k)/dtime
1910	ENDDO
1911	ENDDO
1912	IF (check) THEN
1913	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1914	WRITE(lunout,*) "avantcon=", za
1915	ENDIF
1916	zx_ajustq = .FALSE.
1917	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1918	IF (zx_ajustq) THEN
1919	DO i = 1, klon
1920	z_avant(i) = 0.0
1921	ENDDO
1922	DO k = 1, klev
1923	DO i = 1, klon
1924	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1925	. *(paprs(i,k)-paprs(i,k+1))/RG
1926	ENDDO
1927	ENDDO
1928	ENDIF
1929
1930	c Calcule de vitesse verticale a partir de flux de masse verticale
1931	DO k = 1, klev
1932	DO i = 1, klon
1933	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1934	END DO
1935	END DO
1936
1937	IF (iflag_con.EQ.1) THEN
1938	stop'reactiver le call conlmd dans physiq.F'
1939	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1940	c . d_t_con, d_q_con,
1941	c . rain_con, snow_con, ibas_con, itop_con)
1942	ELSE IF (iflag_con.EQ.2) THEN
1943	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1944	e conv_t, conv_q, -evap, omega,
1945	s d_t_con, d_q_con, rain_con, snow_con,
1946	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1947	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1948	d_u_con = 0.
1949	d_v_con = 0.
1950
1951	WHERE (rain_con < 0.) rain_con = 0.
1952	WHERE (snow_con < 0.) snow_con = 0.
1953	DO i = 1, klon
1954	ibas_con(i) = klev+1 - kcbot(i)
1955	itop_con(i) = klev+1 - kctop(i)
1956	ENDDO
1957	ELSE IF (iflag_con.GE.3) THEN
1958	c nb of tracers for the KE convection:
1959	c MAF la partie traceurs est faite dans phytrac
1960	c on met ntra=1 pour limiter les appels mais on peut
1961	c supprimer les calculs / ftra.
1962	ntra = 1
1963
1964	c=====================================================================================
1965	cajout pour la parametrisation des poches froides:
1966	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1967	do k=1,klev
1968	do i=1,klon
1969	if (iflag_wake.eq.1) then
1970	t_wake(i,k) = t_seri(i,k)
1971	. +(1-wake_s(i))*wake_deltat(i,k)
1972	q_wake(i,k) = q_seri(i,k)
1973	. +(1-wake_s(i))*wake_deltaq(i,k)
1974	t_undi(i,k) = t_seri(i,k)
1975	. -wake_s(i)*wake_deltat(i,k)
1976	q_undi(i,k) = q_seri(i,k)
1977	. -wake_s(i)*wake_deltaq(i,k)
1978	else
1979	t_wake(i,k) = t_seri(i,k)
1980	q_wake(i,k) = q_seri(i,k)
1981	t_undi(i,k) = t_seri(i,k)
1982	q_undi(i,k) = q_seri(i,k)
1983	endif
1984	enddo
1985	enddo
1986
1987	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
1988	cc-- pour le soulevement des particules dans le modele convectif
1989	c
1990	do i = 1,klon
1991	ALE(i) = 0.
1992	ALP(i) = 0.
1993	enddo
1994	c
1995	ccalcul de ale_wake et alp_wake
1996	do i = 1,klon
1997	if (iflag_wake.eq.1) then
1998	ale_wake(i) = 0.5wake_cstar(i)*2
1999	alp_wake(i) = wake_fip(i)
2000	else
2001	ale_wake(i) = 0.
2002	alp_wake(i) = 0.
2003	endif
2004	enddo
2005	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2006	cdans le thermique sinon
2007	if (iflag_coupl.eq.0) then
2008	if (debut) print*,'ALE et ALP imposes'
2009	do i = 1,klon
2010	con ne couple que ale
2011	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2012	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2013	con ne couple que alp
2014	c ALP(i) = alp_wake(i) + Alp_bl(i)
2015	ALP(i) = alp_wake(i) + alp_bl_prescr
2016	enddo
2017	else
2018	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2019	do i = 1,klon
2020	ALE(i) = max(ale_wake(i),Ale_bl(i))
2021	ALP(i) = alp_wake(i) + Alp_bl(i)
2022	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2023	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2024	enddo
2025	endif
2026	do i=1,klon
2027	if (alp(i)>alp_max) then
2028	print*,'WARNING SUPER ALP (seuil=',alp_max,
2029	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2030	alp(i)=alp_max
2031	endif
2032	if (ale(i)>ale_max) then
2033	print*,'WARNING SUPER ALE (seuil=',ale_max,
2034	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2035	ale(i)=ale_max
2036	endif
2037	enddo
2038
2039	cfin calcul ale et alp
2040	c=================================================================================================
2041
2042
2043	c sb, oct02:
2044	c Schema de convection modularise et vectorise:
2045	c (driver commun aux versions 3 et 4)
2046	c
2047	IF (ok_cvl) THEN ! new driver for convectL
2048
2049	CALL concvl (iflag_con,iflag_clos,
2050	. dtime,paprs,pplay,t_undi,q_undi,
2051	. t_wake,q_wake,
2052	. u_seri,v_seri,tr_seri,nbtr,
2053	. ALE,ALP,
2054	. ema_work1,ema_work2,
2055	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2056	. rain_con, snow_con, ibas_con, itop_con, sigd,
2057	. upwd,dnwd,dnwd0,
2058	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2059	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2060	. pmflxr,pmflxs,da,phi,mp,
2061	. ftd,fqd,lalim_conv,wght_th)
2062
2063	cIM begin
2064	print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2065	.dnwd0(1,1),ftd(1,1),fqd(1,1)
2066	cIM end
2067	cIM cf. FH
2068	clwcon0=qcondc
2069	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2070
2071	ELSE ! ok_cvl
2072	c MAF conema3 ne contient pas les traceurs
2073	CALL conema3 (dtime,
2074	. paprs,pplay,t_seri,q_seri,
2075	. u_seri,v_seri,tr_seri,ntra,
2076	. ema_work1,ema_work2,
2077	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2078	. rain_con, snow_con, ibas_con, itop_con,
2079	. upwd,dnwd,dnwd0,bas,top,
2080	. Ma,cape,tvp,rflag,
2081	. pbase
2082	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2083	. ,clwcon0)
2084
2085	ENDIF ! ok_cvl
2086
2087	c
2088	c Correction precip
2089	rain_con = rain_con * cvl_corr
2090	snow_con = snow_con * cvl_corr
2091	c
2092
2093	IF (.NOT. ok_gust) THEN
2094	do i = 1, klon
2095	wd(i)=0.0
2096	enddo
2097	ENDIF
2098
2099	c =================================================================== c
2100	c Calcul des proprietes des nuages convectifs
2101	c
2102
2103	c calcul des proprietes des nuages convectifs
2104	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2105	call clouds_gno
2106	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2107
2108	c =================================================================== c
2109
2110	DO i = 1, klon
2111	ema_pcb(i) = pbase(i)
2112	ENDDO
2113	DO i = 1, klon
2114
2115	! L'idicage de itop_con peut cacher un pb potentiel
2116	! FH sous la dictee de JYG, CR
2117	ema_pct(i) = paprs(i,itop_con(i)+1)
2118
2119	if (itop_con(i).gt.klev-3) then
2120	print*,'La convection monte trop haut '
2121	print*,'itop_con(,',i,',)=',itop_con(i)
2122	endif
2123	ENDDO
2124	DO i = 1, klon
2125	ema_cbmf(i) = ema_workcbmf(i)
2126	ENDDO
2127	ELSE IF (iflag_con.eq.0) THEN
2128	write(lunout,*) 'On n appelle pas la convection'
2129	clwcon0=0.
2130	rnebcon0=0.
2131	d_t_con=0.
2132	d_q_con=0.
2133	d_u_con=0.
2134	d_v_con=0.
2135	rain_con=0.
2136	snow_con=0.
2137	bas=1
2138	top=1
2139	ELSE
2140	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2141	CALL abort
2142	ENDIF
2143
2144	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2145	c . d_u_con, d_v_con)
2146
2147	!-----------------------------------------------------------------------------------------
2148	! ajout des tendances de la diffusion turbulente
2149	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2150	!-----------------------------------------------------------------------------------------
2151
2152	if (mydebug) then
2153	call writefield_phy('u_seri',u_seri,llm)
2154	call writefield_phy('v_seri',v_seri,llm)
2155	call writefield_phy('t_seri',t_seri,llm)
2156	call writefield_phy('q_seri',q_seri,llm)
2157	endif
2158
2159	cIM
2160	IF (ip_ebil_phy.ge.2) THEN
2161	ztit='after convect'
2162	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2163	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2164	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2165	call diagphy(airephy,ztit,ip_ebil_phy
2166	e , zero_v, zero_v, zero_v, zero_v, zero_v
2167	e , zero_v, rain_con, snow_con, ztsol
2168	e , d_h_vcol, d_qt, d_ec
2169	s , fs_bound, fq_bound )
2170	END IF
2171	C
2172	IF (check) THEN
2173	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2174	WRITE(lunout,*)"aprescon=", za
2175	zx_t = 0.0
2176	za = 0.0
2177	DO i = 1, klon
2178	za = za + airephy(i)/FLOAT(klon)
2179	zx_t = zx_t + (rain_con(i)+
2180	. snow_con(i))*airephy(i)/FLOAT(klon)
2181	ENDDO
2182	zx_t = zx_t/za*dtime
2183	WRITE(lunout,*)"Precip=", zx_t
2184	ENDIF
2185	IF (zx_ajustq) THEN
2186	DO i = 1, klon
2187	z_apres(i) = 0.0
2188	ENDDO
2189	DO k = 1, klev
2190	DO i = 1, klon
2191	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2192	. *(paprs(i,k)-paprs(i,k+1))/RG
2193	ENDDO
2194	ENDDO
2195	DO i = 1, klon
2196	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2197	. /z_apres(i)
2198	ENDDO
2199	DO k = 1, klev
2200	DO i = 1, klon
2201	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2202	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2203	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2204	ENDIF
2205	ENDDO
2206	ENDDO
2207	ENDIF
2208	zx_ajustq=.FALSE.
2209
2210	c
2211	c=============================================================================
2212	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2213	cpour la couche limite diffuse pour l instant
2214	c
2215	if (iflag_wake.eq.1) then
2216	DO k=1,klev
2217	DO i=1,klon
2218	dt_dwn(i,k) = ftd(i,k)
2219	wdt_PBL(i,k) = 0.
2220	dq_dwn(i,k) = fqd(i,k)
2221	wdq_PBL(i,k) = 0.
2222	M_dwn(i,k) = dnwd0(i,k)
2223	M_up(i,k) = upwd(i,k)
2224	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2225	udt_PBL(i,k) = 0.
2226	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2227	udq_PBL(i,k) = 0.
2228	ENDDO
2229	ENDDO
2230	c
2231	ccalcul caracteristiques de la poche froide
2232	call calWAKE (paprs,pplay,dtime
2233	: ,t_seri,q_seri,omega
2234	: ,dt_dwn,dq_dwn,M_dwn,M_up
2235	: ,dt_a,dq_a,sigd
2236	: ,wdt_PBL,wdq_PBL
2237	: ,udt_PBL,udq_PBL
2238	o ,wake_deltat,wake_deltaq,wake_dth
2239	o ,wake_h,wake_s,wake_dens
2240	o ,wake_pe,wake_fip,wake_gfl
2241	o ,dt_wake,dq_wake
2242	o ,wake_k, t_undi,q_undi
2243	o ,wake_omgbdth,wake_dp_omgb
2244	o ,wake_dtKE,wake_dqKE
2245	o ,wake_dtPBL,wake_dqPBL
2246	o ,wake_omg,wake_dp_deltomg
2247	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2248	o ,wake_ddeltat,wake_ddeltaq)
2249	c
2250	!-----------------------------------------------------------------------------------------
2251	! ajout des tendances des poches froides
2252	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2253	! coherent avec les autres d_t_...
2254	d_t_wake(:,:)=dt_wake(:,:)*dtime
2255	d_q_wake(:,:)=dq_wake(:,:)*dtime
2256	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2257	!-----------------------------------------------------------------------------------------
2258
2259	endif
2260	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2261	c
2262	c===================================================================
2263	c Convection seche (thermiques ou ajustement)
2264	c===================================================================
2265	c
2266	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2267	s ,seuil_inversion,weak_inversion,dthmin)
2268
2269
2270
2271	d_t_ajsb(:,:)=0.
2272	d_q_ajsb(:,:)=0.
2273	d_t_ajs(:,:)=0.
2274	d_u_ajs(:,:)=0.
2275	d_v_ajs(:,:)=0.
2276	d_q_ajs(:,:)=0.
2277	clwcon0th(:,:)=0.
2278	c
2279	fm_therm(:,:)=0.
2280	entr_therm(:,:)=0.
2281	detr_therm(:,:)=0.
2282	c
2283	IF(prt_level>9)WRITE(lunout,*)
2284	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2285	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2286	if(iflag_thermals.lt.0) then
2287	c Rien
2288	c ====
2289	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2290
2291
2292	else
2293
2294	c Thermiques
2295	c ==========
2296	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2297	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2298
2299
2300	if (iflag_thermals.gt.1) then
2301	call calltherm(pdtphys
2302	s ,pplay,paprs,pphi,weak_inversion
2303	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2304	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2305	s ,fm_therm,entr_therm,detr_therm
2306	s ,zqasc,clwcon0th,lmax_th,ratqscth
2307	s ,ratqsdiff,zqsatth
2308	con rajoute ale et alp, et les caracteristiques de la couche alim
2309	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0)
2310	endif
2311
2312	! call calltherm(pdtphys
2313	! s ,pplay,paprs,pphi
2314	! s ,u_seri,v_seri,t_seri,q_seri
2315	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2316	! s ,fm_therm,entr_therm)
2317
2318	c Ajustement sec
2319	c ==============
2320
2321	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2322	! a partir du sommet des thermiques.
2323	! Dans le cas contraire, on demarre au niveau 1.
2324
2325	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2326
2327	if(iflag_thermals.eq.0) then
2328	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2329	limbas(:)=1
2330	else
2331	limbas(:)=lmax_th(:)
2332	endif
2333
2334	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2335	! pour des test de convergence numerique.
2336	! Le nouveau ajsec est a priori mieux, meme pour le cas
2337	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2338	! non nulles numeriquement pour des mailles non concernees.
2339
2340	if (iflag_thermals.eq.0) then
2341	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2342	s , d_t_ajsb, d_q_ajsb)
2343	else
2344	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2345	s , d_t_ajsb, d_q_ajsb)
2346	endif
2347
2348	!-----------------------------------------------------------------------------------------
2349	! ajout des tendances de l'ajustement sec ou des thermiques
2350	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2351	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2352	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2353
2354	!-----------------------------------------------------------------------------------------
2355
2356	endif
2357
2358	endif
2359	c
2360	c===================================================================
2361	cIM
2362	IF (ip_ebil_phy.ge.2) THEN
2363	ztit='after dry_adjust'
2364	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2365	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2366	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2367	END IF
2368
2369
2370	c-------------------------------------------------------------------------
2371	c Caclul des ratqs
2372	c-------------------------------------------------------------------------
2373
2374	c print*,'calcul des ratqs'
2375	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2376	c ----------------
2377	c on ecrase le tableau ratqsc calcule par clouds_gno
2378	if (iflag_cldcon.eq.1) then
2379	do k=1,klev
2380	do i=1,klon
2381	if(ptconv(i,k)) then
2382	ratqsc(i,k)=ratqsbas
2383	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2384	else
2385	ratqsc(i,k)=0.
2386	endif
2387	enddo
2388	enddo
2389	else if (iflag_cldcon.eq.4) then
2390	ptconvth(:,:)=.false.
2391	ratqsc(:,:)=0.
2392	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2393	call clouds_gno
2394	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2395	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2396	endif
2397
2398	c ratqs stables
2399	c -------------
2400
2401	if (iflag_ratqs.eq.0) then
2402
2403	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2404	do k=1,klev
2405	do i=1, klon
2406	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2407	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2408	enddo
2409	enddo
2410
2411	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2412	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2413	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2414	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2415	! Il s'agit de differents tests dans la phase de reglage du modele
2416	! avec thermiques.
2417
2418	else if (iflag_ratqs.eq.1) then
2419
2420	do k=1,klev
2421	do i=1, klon
2422	if (pplay(i,k).ge.60000.) then
2423	ratqss(i,k)=ratqsbas
2424	else if ((pplay(i,k).ge.30000.).and.
2425	s (pplay(i,k).lt.60000.)) then
2426	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2427	s (60000.-pplay(i,k))/(60000.-30000.)
2428	else
2429	ratqss(i,k)=ratqshaut
2430	endif
2431	enddo
2432	enddo
2433
2434	else
2435
2436	do k=1,klev
2437	do i=1, klon
2438	if (pplay(i,k).ge.60000.) then
2439	ratqss(i,k)=ratqsbas
2440	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2441	else if ((pplay(i,k).ge.30000.).and.
2442	s (pplay(i,k).lt.60000.)) then
2443	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2444	s (60000.-pplay(i,k))/(60000.-30000.)
2445	else
2446	ratqss(i,k)=ratqshaut
2447	endif
2448	enddo
2449	enddo
2450	endif
2451
2452
2453
2454
2455	c ratqs final
2456	c -----------
2457
2458	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2459	s .or.iflag_cldcon.eq.4) then
2460
2461	! On ajoute une constante au ratqsc*2 pour tenir compte de
2462	! fluctuations turbulentes de petite echelle
2463
2464	do k=1,klev
2465	do i=1,klon
2466	if ((fm_therm(i,k).gt.1.e-10)) then
2467	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2468	endif
2469	enddo
2470	enddo
2471
2472	! les ratqs sont une conbinaison de ratqss et ratqsc
2473	! 1800s, en dur pour le moment, est le temps de
2474	! relaxation des ratqs
2475
2476	facteur=exp(-pdtphys/1800.)
2477
2478	print*,'WARNING ratqs a revoir '
2479	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2480	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2481
2482	else
2483	! on ne prend que le ratqs stable pour fisrtilp
2484	ratqs(:,:)=ratqss(:,:)
2485	endif
2486
2487
2488	c
2489	c Appeler le processus de condensation a grande echelle
2490	c et le processus de precipitation
2491	c-------------------------------------------------------------------------
2492	CALL fisrtilp(dtime,paprs,pplay,
2493	. t_seri, q_seri,ptconv,ratqs,
2494	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2495	. rain_lsc, snow_lsc,
2496	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2497	. frac_impa, frac_nucl,
2498	. prfl, psfl, rhcl)
2499
2500	WHERE (rain_lsc < 0) rain_lsc = 0.
2501	WHERE (snow_lsc < 0) snow_lsc = 0.
2502	!-----------------------------------------------------------------------------------------
2503	! ajout des tendances de la diffusion turbulente
2504	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2505	!-----------------------------------------------------------------------------------------
2506	DO k = 1, klev
2507	DO i = 1, klon
2508	cldfra(i,k) = rneb(i,k)
2509	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2510	ENDDO
2511	ENDDO
2512	IF (check) THEN
2513	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2514	WRITE(lunout,*)"apresilp=", za
2515	zx_t = 0.0
2516	za = 0.0
2517	DO i = 1, klon
2518	za = za + airephy(i)/FLOAT(klon)
2519	zx_t = zx_t + (rain_lsc(i)
2520	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2521	ENDDO
2522	zx_t = zx_t/za*dtime
2523	WRITE(lunout,*)"Precip=", zx_t
2524	ENDIF
2525	cIM
2526	IF (ip_ebil_phy.ge.2) THEN
2527	ztit='after fisrt'
2528	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2529	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2530	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2531	call diagphy(airephy,ztit,ip_ebil_phy
2532	e , zero_v, zero_v, zero_v, zero_v, zero_v
2533	e , zero_v, rain_lsc, snow_lsc, ztsol
2534	e , d_h_vcol, d_qt, d_ec
2535	s , fs_bound, fq_bound )
2536	END IF
2537
2538	if (mydebug) then
2539	call writefield_phy('u_seri',u_seri,llm)
2540	call writefield_phy('v_seri',v_seri,llm)
2541	call writefield_phy('t_seri',t_seri,llm)
2542	call writefield_phy('q_seri',q_seri,llm)
2543	endif
2544
2545	c
2546	c-------------------------------------------------------------------
2547	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2548	c-------------------------------------------------------------------
2549
2550	c 1. NUAGES CONVECTIFS
2551	c
2552	cIM cf FH
2553	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2554	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2555	snow_tiedtke=0.
2556	c print*,'avant calcul de la pseudo precip '
2557	c print*,'iflag_cldcon',iflag_cldcon
2558	if (iflag_cldcon.eq.-1) then
2559	rain_tiedtke=rain_con
2560	else
2561	c print*,'calcul de la pseudo precip '
2562	rain_tiedtke=0.
2563	c print*,'calcul de la pseudo precip 0'
2564	do k=1,klev
2565	do i=1,klon
2566	if (d_q_con(i,k).lt.0.) then
2567	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2568	s *(paprs(i,k)-paprs(i,k+1))/rg
2569	endif
2570	enddo
2571	enddo
2572	endif
2573	c
2574	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2575	c
2576
2577	c Nuages diagnostiques pour Tiedtke
2578	CALL diagcld1(paprs,pplay,
2579	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2580	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2581	. diafra,dialiq)
2582	DO k = 1, klev
2583	DO i = 1, klon
2584	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2585	cldliq(i,k) = dialiq(i,k)
2586	cldfra(i,k) = diafra(i,k)
2587	ENDIF
2588	ENDDO
2589	ENDDO
2590
2591	ELSE IF (iflag_cldcon.ge.3) THEN
2592	c On prend pour les nuages convectifs le max du calcul de la
2593	c convection et du calcul du pas de temps precedent diminue d'un facteur
2594	c facttemps
2595	facteur = pdtphys *facttemps
2596	do k=1,klev
2597	do i=1,klon
2598	rnebcon(i,k)=rnebcon(i,k)*facteur
2599	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2600	s then
2601	rnebcon(i,k)=rnebcon0(i,k)
2602	clwcon(i,k)=clwcon0(i,k)
2603	endif
2604	enddo
2605	enddo
2606
2607	c
2608	cjq - introduce the aerosol direct and first indirect radiative forcings
2609	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2610	IF (ok_ade.OR.ok_aie) THEN
2611	IF ( .NOT. aerosol_couple ) THEN
2612	! Get sulfate aerosol distribution
2613	CALL readsulfate(rjourvrai, debut, sulfate)
2614	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2615
2616	! Calculate aerosol optical properties (Olivier Boucher)
2617	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2618	. tau_ae, piz_ae, cg_ae, aerindex)
2619	ENDIF
2620	ELSE
2621	tau_ae(:,:,:)=0.0
2622	piz_ae(:,:,:)=0.0
2623	cg_ae(:,:,:)=0.0
2624	ENDIF
2625
2626	c
2627	cIM calcul nuages par le simulateur ISCCP
2628	c
2629	#ifdef histISCCP
2630	IF (ok_isccp) THEN
2631	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2632	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2633	#include "calcul_simulISCCP.h"
2634	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2635	ENDIF !ok_isccp
2636	#endif
2637
2638	c On prend la somme des fractions nuageuses et des contenus en eau
2639	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2640	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2641
2642	ENDIF
2643
2644	c
2645	c 2. NUAGES STARTIFORMES
2646	c
2647	IF (ok_stratus) THEN
2648	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2649	DO k = 1, klev
2650	DO i = 1, klon
2651	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2652	cldliq(i,k) = dialiq(i,k)
2653	cldfra(i,k) = diafra(i,k)
2654	ENDIF
2655	ENDDO
2656	ENDDO
2657	ENDIF
2658	c
2659	c Precipitation totale
2660	c
2661	DO i = 1, klon
2662	rain_fall(i) = rain_con(i) + rain_lsc(i)
2663	snow_fall(i) = snow_con(i) + snow_lsc(i)
2664	ENDDO
2665	cIM
2666	IF (ip_ebil_phy.ge.2) THEN
2667	ztit="after diagcld"
2668	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2669	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2670	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2671	END IF
2672	c
2673	c Calculer l'humidite relative pour diagnostique
2674	c
2675	DO k = 1, klev
2676	DO i = 1, klon
2677	zx_t = t_seri(i,k)
2678	IF (thermcep) THEN
2679	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2680	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2681	zx_qs = MIN(0.5,zx_qs)
2682	zcor = 1./(1.-retv*zx_qs)
2683	zx_qs = zx_qs*zcor
2684	ELSE
2685	IF (zx_t.LT.t_coup) THEN
2686	zx_qs = qsats(zx_t)/pplay(i,k)
2687	ELSE
2688	zx_qs = qsatl(zx_t)/pplay(i,k)
2689	ENDIF
2690	ENDIF
2691	zx_rh(i,k) = q_seri(i,k)/zx_qs
2692	zqsat(i,k)=zx_qs
2693	ENDDO
2694	ENDDO
2695
2696	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2697	c equivalente a 2m (tpote) pour diagnostique
2698	c
2699	DO i = 1, klon
2700	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2701	IF (thermcep) THEN
2702	IF(zt2m(i).LT.RTT) then
2703	Lheat=RLSTT
2704	ELSE
2705	Lheat=RLVTT
2706	ENDIF
2707	ELSE
2708	IF (zt2m(i).LT.RTT) THEN
2709	Lheat=RLSTT
2710	ELSE
2711	Lheat=RLVTT
2712	ENDIF
2713	ENDIF
2714	tpote(i) = tpot(i)*
2715	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2716	ENDDO
2717
2718	IF (config_inca /= 'none') THEN
2719	#ifdef INCA
2720	CALL VTe(VTphysiq)
2721	CALL VTb(VTinca)
2722	calday = FLOAT(julien) + gmtime
2723
2724	IF (config_inca == 'aero') THEN
2725	CALL AEROSOL_METEO_CALC(
2726	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2727	$ prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2728	END IF
2729
2730	zxsnow_dummy(:) = 0.0
2731
2732	CALL chemhook_begin (calday,
2733	$ julien,
2734	$ gmtime,
2735	$ pctsrf(1,1),
2736	$ rlat,
2737	$ rlon,
2738	$ airephy,
2739	$ paprs,
2740	$ pplay,
2741	$ ycoefh,
2742	$ pphi,
2743	$ t_seri,
2744	$ u,
2745	$ v,
2746	$ wo,
2747	$ q_seri,
2748	$ zxtsol,
2749	$ zxsnow_dummy,
2750	$ solsw,
2751	$ albsol1,
2752	$ rain_fall,
2753	$ snow_fall,
2754	$ itop_con,
2755	$ ibas_con,
2756	$ cldfra,
2757	$ iim,
2758	$ jjm,
2759	$ tr_seri,
2760	$ ftsol,
2761	$ paprs,
2762	$ cdragh,
2763	$ cdragm,
2764	$ pctsrf,
2765	$ pdtphys,
2766	$ itap)
2767
2768	CALL VTe(VTinca)
2769	CALL VTb(VTphysiq)
2770	#endif
2771	END IF !config_inca /= 'none'
2772	c
2773	c Calculer les parametres optiques des nuages et quelques
2774	c parametres pour diagnostiques:
2775	c
2776
2777	IF (aerosol_couple) THEN
2778	sulfate(:,:) = ccm(:,:,1)
2779	sulfate_pi(:,:) = ccm(:,:,2)
2780	ENDIF
2781
2782	if (ok_newmicro) then
2783	CALL newmicro (paprs, pplay,ok_newmicro,
2784	. t_seri, cldliq, cldfra, cldtau, cldemi,
2785	. cldh, cldl, cldm, cldt, cldq,
2786	. flwp, fiwp, flwc, fiwc,
2787	e ok_aie,
2788	e sulfate, sulfate_pi,
2789	e bl95_b0, bl95_b1,
2790	s cldtaupi, re, fl)
2791	else
2792	CALL nuage (paprs, pplay,
2793	. t_seri, cldliq, cldfra, cldtau, cldemi,
2794	. cldh, cldl, cldm, cldt, cldq,
2795	e ok_aie,
2796	e sulfate, sulfate_pi,
2797	e bl95_b0, bl95_b1,
2798	s cldtaupi, re, fl)
2799
2800	endif
2801	c
2802	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2803	c
2804	IF (MOD(itaprad,radpas).EQ.0) THEN
2805
2806	DO i = 1, klon
2807	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2808	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2809	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2810	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2811	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2812	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2813	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2814	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2815	ENDDO
2816
2817	if (mydebug) then
2818	call writefield_phy('u_seri',u_seri,llm)
2819	call writefield_phy('v_seri',v_seri,llm)
2820	call writefield_phy('t_seri',t_seri,llm)
2821	call writefield_phy('q_seri',q_seri,llm)
2822	endif
2823
2824	IF (aerosol_couple) THEN
2825	#ifdef INCA
2826	CALL radlwsw_inca
2827	e (kdlon,kflev,dist, rmu0, fract, solaire,
2828	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2829	e wo,
2830	e cldfra, cldemi, cldtau,
2831	s heat,heat0,cool,cool0,radsol,albpla,
2832	s topsw,toplw,solsw,sollw,
2833	s sollwdown,
2834	s topsw0,toplw0,solsw0,sollw0,
2835	s lwdn0, lwdn, lwup0, lwup,
2836	s swdn0, swdn, swup0, swup,
2837	e ok_ade, ok_aie,
2838	e tau_inca, piz_inca, cg_inca,
2839	s topswad_inca, solswad_inca,
2840	s topswad0_inca, solswad0_inca,
2841	s topsw_inca, topsw0_inca,
2842	s solsw_inca, solsw0_inca,
2843	e cldtaupi,
2844	s topswai_inca, solswai_inca)
2845	#endif
2846	ELSE
2847	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2848	e (dist, rmu0, fract,
2849	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2850	e wo,
2851	e cldfra, cldemi, cldtau,
2852	s heat,heat0,cool,cool0,radsol,albpla,
2853	s topsw,toplw,solsw,sollw,
2854	s sollwdown,
2855	s topsw0,toplw0,solsw0,sollw0,
2856	s lwdn0, lwdn, lwup0, lwup,
2857	s swdn0, swdn, swup0, swup,
2858	e ok_ade, ok_aie, ! new for aerosol radiative effects
2859	e tau_ae, piz_ae, cg_ae, ! ="=
2860	s topswad, solswad, ! ="=
2861	e cldtaupi, ! ="=
2862	s topswai, solswai,zqsat,flwc,fiwc) ! ="=
2863	ENDIF
2864	itaprad = 0
2865	ENDIF
2866	itaprad = itaprad + 1
2867
2868	if (iflag_radia.eq.0) then
2869	print *,'--------------------------------------------------'
2870	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2871	print *,'>>>> heat et cool mis a zero '
2872	print *,'--------------------------------------------------'
2873	heat=0.
2874	cool=0.
2875	endif
2876
2877
2878	c
2879	c Ajouter la tendance des rayonnements (tous les pas)
2880	c
2881	DO k = 1, klev
2882	DO i = 1, klon
2883	t_seri(i,k) = t_seri(i,k)
2884	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2885	ENDDO
2886	ENDDO
2887	c
2888	if (mydebug) then
2889	call writefield_phy('u_seri',u_seri,llm)
2890	call writefield_phy('v_seri',v_seri,llm)
2891	call writefield_phy('t_seri',t_seri,llm)
2892	call writefield_phy('q_seri',q_seri,llm)
2893	endif
2894
2895	cIM
2896	IF (ip_ebil_phy.ge.2) THEN
2897	ztit='after rad'
2898	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2899	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2900	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2901	call diagphy(airephy,ztit,ip_ebil_phy
2902	e , topsw, toplw, solsw, sollw, zero_v
2903	e , zero_v, zero_v, zero_v, ztsol
2904	e , d_h_vcol, d_qt, d_ec
2905	s , fs_bound, fq_bound )
2906	END IF
2907	c
2908	c
2909	c Calculer l'hydrologie de la surface
2910	c
2911	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2912	c . agesno, ftsol,fqsurf,fsnow, ruis)
2913	c
2914
2915	c
2916	c Calculer le bilan du sol et la derive de temperature (couplage)
2917	c
2918	DO i = 1, klon
2919	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2920	c a la demande de JLD
2921	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2922	ENDDO
2923	c
2924	cmoddeblott(jan95)
2925	c Appeler le programme de parametrisation de l'orographie
2926	c a l'echelle sous-maille:
2927	c
2928	IF (ok_orodr) THEN
2929	c
2930	c selection des points pour lesquels le shema est actif:
2931	igwd=0
2932	DO i=1,klon
2933	itest(i)=0
2934	c IF ((zstd(i).gt.10.0)) THEN
2935	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2936	itest(i)=1
2937	igwd=igwd+1
2938	idx(igwd)=i
2939	ENDIF
2940	ENDDO
2941	c igwdim=MAX(1,igwd)
2942	c
2943	IF (ok_strato) THEN
2944
2945	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
2946	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2947	e igwd,idx,itest,
2948	e t_seri, u_seri, v_seri,
2949	s zulow, zvlow, zustrdr, zvstrdr,
2950	s d_t_oro, d_u_oro, d_v_oro)
2951
2952	ELSE
2953	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2954	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2955	e igwd,idx,itest,
2956	e t_seri, u_seri, v_seri,
2957	s zulow, zvlow, zustrdr, zvstrdr,
2958	s d_t_oro, d_u_oro, d_v_oro)
2959	ENDIF
2960	c
2961	c ajout des tendances
2962	!-----------------------------------------------------------------------------------------
2963	! ajout des tendances de la trainee de l'orographie
2964	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
2965	!-----------------------------------------------------------------------------------------
2966	c
2967	ENDIF ! fin de test sur ok_orodr
2968	c
2969	if (mydebug) then
2970	call writefield_phy('u_seri',u_seri,llm)
2971	call writefield_phy('v_seri',v_seri,llm)
2972	call writefield_phy('t_seri',t_seri,llm)
2973	call writefield_phy('q_seri',q_seri,llm)
2974	endif
2975
2976	IF (ok_orolf) THEN
2977	c
2978	c selection des points pour lesquels le shema est actif:
2979	igwd=0
2980	DO i=1,klon
2981	itest(i)=0
2982	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2983	itest(i)=1
2984	igwd=igwd+1
2985	idx(igwd)=i
2986	ENDIF
2987	ENDDO
2988	c igwdim=MAX(1,igwd)
2989	c
2990	IF (ok_strato) THEN
2991
2992	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
2993	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
2994	e igwd,idx,itest,
2995	e t_seri, u_seri, v_seri,
2996	s zulow, zvlow, zustrli, zvstrli,
2997	s d_t_lif, d_u_lif, d_v_lif )
2998
2999	ELSE
3000	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3001	e rlat,zmea,zstd,zpic,
3002	e itest,
3003	e t_seri, u_seri, v_seri,
3004	s zulow, zvlow, zustrli, zvstrli,
3005	s d_t_lif, d_u_lif, d_v_lif)
3006	ENDIF
3007	c
3008	!-----------------------------------------------------------------------------------------
3009	! ajout des tendances de la portance de l'orographie
3010	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3011	!-----------------------------------------------------------------------------------------
3012	c
3013	ENDIF ! fin de test sur ok_orolf
3014	C HINES GWD PARAMETRIZATION
3015
3016	IF (ok_hines) then
3017
3018	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3019	i rlat,t_seri,u_seri,v_seri,
3020	o zustrhi,zvstrhi,
3021	o d_t_hin, d_u_hin, d_v_hin)
3022	c
3023	c ajout des tendances
3024	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3025
3026	ENDIF
3027	c
3028
3029	c
3030	cIM cf. FLott BEG
3031	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3032
3033	if (mydebug) then
3034	call writefield_phy('u_seri',u_seri,llm)
3035	call writefield_phy('v_seri',v_seri,llm)
3036	call writefield_phy('t_seri',t_seri,llm)
3037	call writefield_phy('q_seri',q_seri,llm)
3038	endif
3039
3040	DO i = 1, klon
3041	zustrph(i)=0.
3042	zvstrph(i)=0.
3043	ENDDO
3044	DO k = 1, klev
3045	DO i = 1, klon
3046	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3047	c (paprs(i,k)-paprs(i,k+1))/rg
3048	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3049	c (paprs(i,k)-paprs(i,k+1))/rg
3050	ENDDO
3051	ENDDO
3052	c
3053	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3054	c
3055	IF (is_sequential) THEN
3056
3057	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3058	C ra,rg,romega,
3059	C rlat,rlon,pphis,
3060	C zustrdr,zustrli,zustrph,
3061	C zvstrdr,zvstrli,zvstrph,
3062	C paprs,u,v,
3063	C aam, torsfc)
3064	ENDIF
3065	cIM cf. FLott END
3066	cIM
3067	IF (ip_ebil_phy.ge.2) THEN
3068	ztit='after orography'
3069	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3070	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3071	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3072	END IF
3073	c
3074	c
3075	cAA
3076	cAA Installation de l'interface online-offline pour traceurs
3077	cAA
3078	c====================================================================
3079	c Calcul des tendances traceurs
3080	c====================================================================
3081	C
3082	IF (config_inca /= 'none') rnpb=.FALSE.
3083
3084	call phytrac ( rnpb,
3085	I itap,
3086	I julien,
3087	I gmtime,
3088	I debut,
3089	I lafin,
3090	I nqmax-2,
3091	I nlon,
3092	I nlev,
3093	I dtime,
3094	I u,
3095	I v,
3096	I t,
3097	I paprs,
3098	I pplay,
3099	I pmfu,
3100	I pmfd,
3101	I pen_u,
3102	I pde_u,
3103	I pen_d,
3104	I pde_d,
3105	I ycoefh,
3106	I fm_therm,
3107	I entr_therm,
3108	I yu1,
3109	I yv1,
3110	I ftsol,
3111	I pctsrf,
3112	I rlat,
3113	I frac_impa,
3114	I frac_nucl,
3115	I rlon,
3116	I presnivs,
3117	I pphis,
3118	I pphi,
3119	I albsol1,
3120	I qx(1,1,1),
3121	I rhcl,
3122	I cldfra,
3123	I rneb,
3124	I diafra,
3125	I cldliq,
3126	I itop_con,
3127	I ibas_con,
3128	I pmflxr,
3129	I pmflxs,
3130	I prfl,
3131	I psfl,
3132	I da,
3133	I phi,
3134	I mp,
3135	I upwd,
3136	I dnwd,
3137	I aerosol_couple,
3138	I flxmass_w,
3139	I tau_inca,
3140	I piz_inca,
3141	I cg_inca,
3142	I ccm,
3143	I rfname,
3144	O tr_seri)
3145
3146	IF (offline) THEN
3147
3148	print*,'Attention on met a 0 les thermiques pour phystoke'
3149	call phystokenc (
3150	I nlon,nlev,pdtphys,rlon,rlat,
3151	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3152	I fm_therm,entr_therm,
3153	I ycoefh,yu1,yv1,ftsol,pctsrf,
3154	I frac_impa, frac_nucl,
3155	I pphis,airephy,dtime,itap)
3156
3157
3158	ENDIF
3159
3160	c
3161	c Calculer le transport de l'eau et de l'energie (diagnostique)
3162	c
3163	CALL transp (paprs,zxtsol,
3164	e t_seri, q_seri, u_seri, v_seri, zphi,
3165	s ve, vq, ue, uq)
3166	c
3167	cIM global posePB BEG
3168	IF(1.EQ.0) THEN
3169	c
3170	CALL transp_lay (paprs,zxtsol,
3171	e t_seri, q_seri, u_seri, v_seri, zphi,
3172	s ve_lay, vq_lay, ue_lay, uq_lay)
3173	c
3174	ENDIF !(1.EQ.0) THEN
3175	cIM global posePB END
3176	c Accumuler les variables a stocker dans les fichiers histoire:
3177	c
3178	c+jld ec_conser
3179	DO k = 1, klev
3180	DO i = 1, klon
3181	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3182	d_t_ec(i,k)=0.5/ZRCPD
3183	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3184	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3185	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3186	END DO
3187	END DO
3188	c-jld ec_conser
3189	cIM
3190	IF (ip_ebil_phy.ge.1) THEN
3191	ztit='after physic'
3192	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3193	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3194	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3195	C Comme les tendances de la physique sont ajoute dans la dynamique,
3196	C on devrait avoir que la variation d'entalpie par la dynamique
3197	C est egale a la variation de la physique au pas de temps precedent.
3198	C Donc la somme de ces 2 variations devrait etre nulle.
3199	call diagphy(airephy,ztit,ip_ebil_phy
3200	e , topsw, toplw, solsw, sollw, sens
3201	e , evap, rain_fall, snow_fall, ztsol
3202	e , d_h_vcol, d_qt, d_ec
3203	s , fs_bound, fq_bound )
3204	C
3205	d_h_vcol_phy=d_h_vcol
3206	C
3207	END IF
3208	C
3209	c=======================================================================
3210	c SORTIES
3211	c=======================================================================
3212
3213	cIM Interpolation sur les niveaux de pression du NMC
3214	c -------------------------------------------------
3215	c
3216	#include "calcul_STDlev.h"
3217	c
3218	c slp sea level pressure
3219	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3220	c
3221	ccc prw = eau precipitable
3222	DO i = 1, klon
3223	prw(i) = 0.
3224	DO k = 1, klev
3225	prw(i) = prw(i) +
3226	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3227	ENDDO
3228	ENDDO
3229	c
3230	cIM initialisation + calculs divers diag AMIP2
3231	c
3232	#include "calcul_divers.h"
3233	c
3234	IF (config_inca /= 'none') THEN
3235	#ifdef INCA
3236	CALL VTe(VTphysiq)
3237	CALL VTb(VTinca)
3238
3239	CALL chemhook_end (calday,
3240	$ dtime,
3241	$ pplay,
3242	$ t_seri,
3243	$ tr_seri,
3244	$ nbtr,
3245	$ paprs,
3246	$ q_seri,
3247	$ annee_ref,
3248	$ day_ini,
3249	$ airephy,
3250	$ xjour,
3251	$ pphi,
3252	$ pphis,
3253	$ zx_rh)
3254
3255	CALL VTe(VTinca)
3256	CALL VTb(VTphysiq)
3257	#endif
3258	END IF
3259
3260	c=============================================================
3261	c
3262	c Convertir les incrementations en tendances
3263	c
3264	if (mydebug) then
3265	call writefield_phy('u_seri',u_seri,llm)
3266	call writefield_phy('v_seri',v_seri,llm)
3267	call writefield_phy('t_seri',t_seri,llm)
3268	call writefield_phy('q_seri',q_seri,llm)
3269	endif
3270
3271	DO k = 1, klev
3272	DO i = 1, klon
3273	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3274	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3275	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3276	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3277	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3278	ENDDO
3279	ENDDO
3280	c
3281	IF (nqmax.GE.3) THEN
3282	DO iq = 3, nqmax
3283	DO k = 1, klev
3284	DO i = 1, klon
3285	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3286	ENDDO
3287	ENDDO
3288	ENDDO
3289	ENDIF
3290	c
3291	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3292	cIM global posePB#include "write_bilKP_ins.h"
3293	cIM global posePB#include "write_bilKP_ave.h"
3294	c
3295	c Sauvegarder les valeurs de t et q a la fin de la physique:
3296	c
3297	DO k = 1, klev
3298	DO i = 1, klon
3299	t_ancien(i,k) = t_seri(i,k)
3300	q_ancien(i,k) = q_seri(i,k)
3301	ENDDO
3302	ENDDO
3303	c
3304	!==========================================================================
3305	! Sorties des tendances pour un point particulier
3306	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3307	! pour le debug
3308	! La valeur de igout est attribuee plus haut dans le programme
3309	!==========================================================================
3310
3311	if (prt_level.ge.1) then
3312	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3313	write(lunout,*)
3314	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
3315	write(lunout,*)
3316	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys,
3317	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3318	s pctsrf(igout,is_sic)
3319	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3320	do k=1,nlev
3321	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3322	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3323	s d_t_eva(igout,k)
3324	enddo
3325	write(lunout,*) 'cool,heat'
3326	do k=1,nlev
3327	write(lunout,*) cool(igout,k),heat(igout,k)
3328	enddo
3329
3330	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3331	do k=1,nlev
3332	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3333	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3334	enddo
3335
3336	write(lunout,*) 'd_ps ',d_ps(igout)
3337	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3338	do k=1,nlev
3339	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3340	s d_qx(igout,k,1),d_qx(igout,k,2)
3341	enddo
3342	endif
3343
3344	!==========================================================================
3345
3346	c============================================================
3347	c Calcul de la temperature potentielle
3348	c============================================================
3349	DO k = 1, klev
3350	DO i = 1, klon
3351	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3352	ENDDO
3353	ENDDO
3354	c
3355
3356	c 22.03.04 BEG
3357	c=============================================================
3358	c Ecriture des sorties
3359	c=============================================================
3360	#ifdef CPP_IOIPSL
3361
3362	c Recupere des varibles calcule dans differents modules
3363	c pour ecriture dans histxxx.nc
3364
3365	! Get some variables from module fonte_neige_mod
3366	CALL fonte_neige_get_vars(pctsrf,
3367	. zxfqcalving, zxfqfonte, zxffonte)
3368
3369
3370	c Commente par abderrahmane le 11 2 08
3371	c#ifdef histhf
3372	c#include "write_histhf.h"
3373	c#endif
3374
3375	c#ifdef histday
3376	c#include "write_histday.h"
3377	c#endif
3378
3379	c#ifdef histmth
3380	c#include "write_histmth.h"
3381	c#endif
3382
3383	c#ifdef histins
3384	c#include "write_histins.h"
3385	c#endif
3386
3387	#include "phys_output_write.h"
3388
3389	#ifdef histISCCP
3390	#include "write_histISCCP.h"
3391	#endif
3392
3393	#ifdef histmthNMC
3394	#include "write_histmthNMC.h"
3395	#endif
3396
3397	#include "write_histday_seri.h"
3398
3399	#include "write_paramLMDZ_phy.h"
3400
3401	#endif
3402
3403	c 22.03.04 END
3404	c
3405	c====================================================================
3406	c Si c'est la fin, il faut conserver l'etat de redemarrage
3407	c====================================================================
3408	c
3409
3410
3411	IF (lafin) THEN
3412	itau_phy = itau_phy + itap
3413	CALL phyredem ("restartphy.nc")
3414	! open(97,form="unformatted",file="finbin")
3415	! write(97) u_seri,v_seri,t_seri,q_seri
3416	! close(97)
3417	ENDIF
3418
3419
3420	RETURN
3421	END
3422	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3423	IMPLICIT none
3424	c
3425	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3426	c la conservation de l'eau
3427	c
3428	#include "YOMCST.h"
3429	INTEGER klon,klev
3430	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3431	REAL aire(klon)
3432	REAL qtotal, zx, qcheck
3433	INTEGER i, k
3434	c
3435	zx = 0.0
3436	DO i = 1, klon
3437	zx = zx + aire(i)
3438	ENDDO
3439	qtotal = 0.0
3440	DO k = 1, klev
3441	DO i = 1, klon
3442	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3443	. *(paprs(i,k)-paprs(i,k+1))/RG
3444	ENDDO
3445	ENDDO
3446	c
3447	qcheck = qtotal/zx
3448	c
3449	RETURN
3450	END
3451	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3452	IMPLICIT none
3453	c
3454	c Tranformer une variable de la grille physique a
3455	c la grille d'ecriture
3456	c
3457	INTEGER nfield,nlon,iim,jjmp1, jjm
3458	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3459	c
3460	INTEGER i, n, ig
3461	c
3462	jjm = jjmp1 - 1
3463	DO n = 1, nfield
3464	DO i=1,iim
3465	ecrit(i,n) = fi(1,n)
3466	ecrit(i+jjm*iim,n) = fi(nlon,n)
3467	ENDDO
3468	DO ig = 1, nlon - 2
3469	ecrit(iim+ig,n) = fi(1+ig,n)
3470	ENDDO
3471	ENDDO
3472	RETURN
3473	END

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

Download in other formats:

Original Format