Context Navigation

physiq.F @ 1037

Last change on this file since 1037 was 1037, checked in by lmdzadmin, 16 years ago
Correction test coherence flags Tiedtke FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 113.3 KB

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

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

Download in other formats:

Original Format