Context Navigation

physiq.F @ 987

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

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

Download in other formats:

Original Format