Context Navigation

physiq.F @ 1054

Last change on this file since 1054 was 1054, checked in by lmdzadmin, 16 years ago
Ajout sorties tendances dynamiques histLES ACA/FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 114.4 KB

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

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

Download in other formats:

Original Format