Context Navigation

physiq.F @ 1010

Last change on this file since 1010 was 1001, checked in by Laurent Fairhead, 16 years ago

Modifs sur le parallelisme: masquage dans la physique
Inclusion strato
mise en coherence etat0
le mode offline fonctionne maintenant en parallele,
les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF

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

File size: 111.0 KB

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

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

Download in other formats:

Original Format