Context Navigation

physiq.F @ 1102

Last change on this file since 1102 was 1102, checked in by idelkadi, 15 years ago
Fichiers run.def gcm.def physiq.def adaptes pour la version actuelle Nettoyage de physiq.F des appels aux ini_hists et write_hists
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.0 KB

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

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

Download in other formats:

Original Format