Context Navigation

physiq.F @ 959

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

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

Download in other formats:

Original Format