Context Navigation

physiq.F @ 1043

Last change on this file since 1043 was 1043, checked in by lmdzadmin, 16 years ago
Corrections sorties KE sous flag iflag_con + choix Tiedtke / sorties ISCCP - incompatibles IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 113.4 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c#define IO_DEBUG
6
7	SUBROUTINE physiq (nlon,nlev,nqmax,
8	. debut,lafin,rjourvrai,gmtime,pdtphys,
9	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
10	. u,v,t,qx,
11	. flxmass_w,
12	. d_u, d_v, d_t, d_qx, d_ps
13	. , dudyn
14	. , PVteta)
15
16	USE ioipsl
17	USE comgeomphy
18	USE write_field_phy
19	USE dimphy
20	USE mod_grid_phy_lmdz
21	USE mod_phys_lmdz_para
22	USE iophy
23	USE misc_mod, mydebug=>debug
24	USE vampir
25	USE pbl_surface_mod, ONLY : pbl_surface
26	USE change_srf_frac_mod
27	USE surface_data, ONLY : type_ocean, ok_veget
28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
29	USE phys_state_var_mod ! Variables sauvegardees de la physique
30	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
31	USE phys_output_mod
32
33	IMPLICIT none
34	c======================================================================
35	c
36	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
37	c
38	c Objet: Moniteur general de la physique du modele
39	cAA Modifications quant aux traceurs :
40	cAA - uniformisation des parametrisations ds phytrac
41	cAA - stockage des moyennes des champs necessaires
42	cAA en mode traceur off-line
43	c======================================================================
44	c CLEFS CPP POUR LES IO
45	c =====================
46	c#define histhf
47	#define histday
48	#define histmth
49	c#define histmthNMC
50	c#define histins
51	c#define histISCCP
52	c======================================================================
53	c modif ( P. Le Van , 12/10/98 )
54	c
55	c Arguments:
56	c
57	c nlon----input-I-nombre de points horizontaux
58	c nlev----input-I-nombre de couches verticales
59	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
60	c debut---input-L-variable logique indiquant le premier passage
61	c lafin---input-L-variable logique indiquant le dernier passage
62	c rjour---input-R-numero du jour de l'experience
63	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
64	c pdtphys-input-R-pas d'integration pour la physique (seconde)
65	c paprs---input-R-pression pour chaque inter-couche (en Pa)
66	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
67	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
68	c pphis---input-R-geopotentiel du sol
69	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
70	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
71	c v-------input-R-vitesse Y (de S a N) en m/s
72	c t-------input-R-temperature (K)
73	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
74	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
75	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
76	c flxmass_w -input-R- flux de masse verticale
77	c d_u-----output-R-tendance physique de "u" (m/s/s)
78	c d_v-----output-R-tendance physique de "v" (m/s/s)
79	c d_t-----output-R-tendance physique de "t" (K/s)
80	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
81	c d_ps----output-R-tendance physique de la pression au sol
82	cIM
83	c PVteta--output-R-vorticite potentielle a des thetas constantes
84	c======================================================================
85	#include "dimensions.h"
86	integer jjmp1
87	parameter (jjmp1=jjm+1-1/jjm)
88	integer iip1
89	parameter (iip1=iim+1)
90
91	#include "regdim.h"
92	#include "indicesol.h"
93	#include "dimsoil.h"
94	#include "clesphys.h"
95	#include "control.h"
96	!#include "logic.h"
97	#include "temps.h"
98	cym#include "comgeomphy.h"
99	#include "advtrac.h"
100	#include "iniprint.h"
101	#include "thermcell.h"
102	c======================================================================
103	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
104	PARAMETER (ok_cvl=.TRUE.)
105	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
106	PARAMETER (ok_gust=.FALSE.)
107	integer iflag_radia ! active ou non le rayonnement (MPL)
108	save iflag_radia
109	c$OMP THREADPRIVATE(iflag_radia)
110	c======================================================================
111	LOGICAL check ! Verifier la conservation du modele en eau
112	PARAMETER (check=.FALSE.)
113	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
114	PARAMETER (ok_stratus=.FALSE.)
115	c======================================================================
116	LOGICAL, SAVE :: rnpb=.TRUE.
117	c$OMP THREADPRIVATE(rnpb)
118	REAL amn, amx
119	INTEGER igout
120	c======================================================================
121	c Clef controlant l'activation du cycle diurne:
122	ccc LOGICAL cycle_diurne
123	ccc PARAMETER (cycle_diurne=.FALSE.)
124	c======================================================================
125	c Modele thermique du sol, a activer pour le cycle diurne:
126	ccc LOGICAL soil_model
127	ccc PARAMETER (soil_model=.FALSE.)
128	c======================================================================
129	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
130	c le calcul du rayonnement est celle apres la precipitation des nuages.
131	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
132	c la condensation et la precipitation. Cette cle augmente les impacts
133	c radiatifs des nuages.
134	ccc LOGICAL new_oliq
135	ccc PARAMETER (new_oliq=.FALSE.)
136	c======================================================================
137	c Clefs controlant deux parametrisations de l'orographie:
138	cc LOGICAL ok_orodr
139	ccc PARAMETER (ok_orodr=.FALSE.)
140	ccc LOGICAL ok_orolf
141	ccc PARAMETER (ok_orolf=.FALSE.)
142	c======================================================================
143	LOGICAL ok_journe ! sortir le fichier journalier
144	save ok_journe
145	c$OMP THREADPRIVATE(ok_journe)
146	c
147	LOGICAL ok_mensuel ! sortir le fichier mensuel
148	save ok_mensuel
149	c$OMP THREADPRIVATE(ok_mensuel)
150	c
151	LOGICAL ok_instan ! sortir le fichier instantane
152	save ok_instan
153	c$OMP THREADPRIVATE(ok_instan)
154	c
155	LOGICAL ok_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 = 'ISCCP-like outputs may be available for KE
1264	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1265	CALL abort_gcm (modname,abort_message,1)
1266	ENDIF
1267	c
1268	c Initialiser les compteurs:
1269	c
1270	itap = 0
1271	itaprad = 0
1272
1273	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1274	!! Un petit travail à faire ici.
1275	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1276
1277	if (iflag_pbl>1) then
1278	PRINT*, "Using method MELLOR&YAMADA"
1279	endif
1280
1281	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1282	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1283	! dyn3d
1284	! Attention : la version precedente n'etait pas tres propre.
1285	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1286	! pour obtenir le meme resultat.
1287	dtime=pdtphys
1288	radpas = NINT( 86400./dtime/nbapp_rad)
1289	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1290
1291	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1292	cIM begin
1293	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1294	$,ratqs(1,1)
1295	cIM end
1296
1297
1298
1299	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1300	c
1301	C on remet le calendrier a zero
1302	c
1303	IF (raz_date .eq. 1) THEN
1304	itau_phy = 0
1305	ENDIF
1306
1307	cIM cf. AM 081204 BEG
1308	PRINT*,'cycle_diurne3 =',cycle_diurne
1309	cIM cf. AM 081204 END
1310	c
1311	CALL printflag( tabcntr0,radpas,ok_journe,
1312	, ok_instan, ok_region )
1313	c
1314	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1315	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1316	. pdtphys
1317	abort_message='Pas physique n est pas correct '
1318	! call abort_gcm(modname,abort_message,1)
1319	dtime=pdtphys
1320	ENDIF
1321	IF (nlon .NE. klon) THEN
1322	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1323	. klon
1324	abort_message='nlon et klon ne sont pas coherents'
1325	call abort_gcm(modname,abort_message,1)
1326	ENDIF
1327	IF (nlev .NE. klev) THEN
1328	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1329	. klev
1330	abort_message='nlev et klev ne sont pas coherents'
1331	call abort_gcm(modname,abort_message,1)
1332	ENDIF
1333	c
1334	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1335	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1336	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1337	abort_message='Nbre d appels au rayonnement insuffisant'
1338	call abort_gcm(modname,abort_message,1)
1339	ENDIF
1340	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1341	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1342	. ok_cvl
1343	c
1344	cKE43
1345	c Initialisation pour la convection de K.E. (sb):
1346	IF (iflag_con.GE.3) THEN
1347
1348	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1349	WRITE(lunout,*)
1350	. "On va utiliser le melange convectif des traceurs qui"
1351	WRITE(lunout,*)"est calcule dans convect4.3"
1352	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1353
1354	DO i = 1, klon
1355	ema_cbmf(i) = 0.
1356	ema_pcb(i) = 0.
1357	ema_pct(i) = 0.
1358	ema_workcbmf(i) = 0.
1359	ENDDO
1360	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1361	DO i = 1, klon
1362	ibas_con(i) = 1
1363	itop_con(i) = 1
1364	ENDDO
1365	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1366	c===============================================================================
1367	cCR:04.12.07: initialisations poches froides
1368	c Controle de ALE et ALP pour la fermeture convective (jyg)
1369	if (iflag_wake.eq.1) then
1370	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1371	s ,alp_bl_prescr, ale_bl_prescr)
1372	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1373	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1374	endif
1375
1376	do i = 1,klon
1377	Ale_bl(i)=0.
1378	Alp_bl(i)=0.
1379	enddo
1380
1381	c================================================================================
1382
1383	ENDIF
1384
1385	rugoro=0.
1386	c34EK
1387	IF (ok_orodr) THEN
1388
1389	rugoro=0.
1390
1391	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1393	! justement quand ok_orodr = false.
1394	! ce rugoro est utilise par la couche limite et fait double emploi
1395	! avec les paramétrisations spécifiques de Francois Lott.
1396	! DO i=1,klon
1397	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1398	! ENDDO
1399	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1400	IF (ok_strato) THEN
1401	CALL SUGWD_strato(klon,klev,paprs,pplay)
1402	ELSE
1403	CALL SUGWD(klon,klev,paprs,pplay)
1404	ENDIF
1405
1406	DO i=1,klon
1407	zuthe(i)=0.
1408	zvthe(i)=0.
1409	if(zstd(i).gt.10.)then
1410	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1411	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1412	endif
1413	ENDDO
1414	ENDIF
1415	c
1416	c
1417	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1418	WRITE(lunout,*)'La frequence de lecture surface est de ',
1419	. lmt_pas
1420	c
1421	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1422	c IF (ok_mensuel) THEN
1423	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1424	c . ecrit_mth
1425	c ENDIF
1426	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1427	c IF (ok_journe) THEN
1428	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1429	c . ecrit_day
1430	c ENDIF
1431	cIM 130904 BEG
1432	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1433	cIM 170305
1434	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1435	cIM 230305
1436	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1437	c
1438	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1439	c
1440	cIM200505 IF (ok_journe) THEN
1441	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1442	cIM200505 . ecrit_hf
1443	cIM200505 ENDIF
1444	cIM 130904 END
1445	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1446	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1447	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1448	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1449	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1450	cIM200505 IF (ok_instan) THEN
1451	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1452	cIM200505 . ecrit_ins
1453	cIM200505 ENDIF
1454	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1455	cIM200505 IF (ok_region) THEN
1456	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1457	cIM200505 . ecrit_reg
1458	cIM200505 ENDIF
1459	cIM 030306 BEG
1460	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1461	cIM : ne pas modifier ecrit_hf2mth
1462	c
1463	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1464	ecrit_hf2mth = ecrit_mth/ecrit_hf
1465	c ecrit_ins en secondes, chaque pas de temps de la physique
1466	ecrit_ins = dtime
1467	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1468	ecrit_hf = ecrit_hf * un_jour
1469	!IM
1470	IF(ecrit_day.LE.1.) THEN
1471	ecrit_day = ecrit_day * un_jour !en secondes
1472	ENDIF
1473	!IM
1474	ecrit_mth = ecrit_mth * un_jour
1475	ecrit_reg = ecrit_reg * un_jour
1476	ecrit_tra = ecrit_tra * un_jour
1477	ecrit_ISCCP = ecrit_ISCCP * un_jour
1478	c
1479	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1480	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1481	. ecrit_hf2mth
1482	cIM 030306 END
1483
1484	capemaxcels = 't_max(X)'
1485	t2mincels = 't_min(X)'
1486	t2maxcels = 't_max(X)'
1487	tinst = 'inst(X)'
1488	tave = 'ave(X)'
1489	cIM cf. AM 081204 BEG
1490	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1491	cIM cf. AM 081204 END
1492	c
1493	c=============================================================
1494	c Initialisation des sorties
1495	c=============================================================
1496
1497	#ifdef CPP_IOIPSL
1498
1499	c Commente par abderrahmane 11 2 08
1500	c#ifdef histhf
1501	c#include "ini_histhf.h"
1502	c#endif
1503
1504	c#ifdef histday
1505	c#include "ini_histday.h"
1506	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1507	c#include "ini_bilKP_ins.h"
1508	c#include "ini_bilKP_ave.h"
1509	c#endif
1510
1511	c#ifdef histmth
1512	c#include "ini_histmth.h"
1513	c#endif
1514
1515	c#ifdef histins
1516	c#include "ini_histins.h"
1517	c#endif
1518
1519	c$OMP MASTER
1520	call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta,
1521	& ctetaSTD,dtime,presnivs,ok_veget,
1522	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1523	& ok_hf,ok_instan,nid_files)
1524	c$OMP END MASTER
1525	c$OMP BARRIER
1526
1527	#ifdef histISCCP
1528	#include "ini_histISCCP.h"
1529	#endif
1530
1531	#ifdef histmthNMC
1532	#include "ini_histmthNMC.h"
1533	#endif
1534
1535	#include "ini_histday_seri.h"
1536
1537	#include "ini_paramLMDZ_phy.h"
1538
1539	#endif
1540
1541	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1542	date0 = zjulian
1543	C date0 = day_ini
1544	WRITE(,) 'physiq date0 : ',date0
1545	c
1546	c
1547	c
1548	c Prescrire l'ozone dans l'atmosphere
1549	c
1550	c
1551	cc DO i = 1, klon
1552	cc DO k = 1, klev
1553	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1554	cc ENDDO
1555	cc ENDDO
1556	c
1557	IF (config_inca /= 'none') THEN
1558	#ifdef INCA
1559	CALL VTe(VTphysiq)
1560	CALL VTb(VTinca)
1561	iii = MOD(NINT(xjour),360)
1562	calday = FLOAT(iii) + gmtime
1563	WRITE(lunout,*) 'initial time ', xjour, calday
1564
1565	CALL chemini(
1566	$ rg,
1567	$ ra,
1568	$ airephy,
1569	$ rlat,
1570	$ rlon,
1571	$ presnivs,
1572	$ calday,
1573	$ klon,
1574	$ nqmax,
1575	$ pdtphys,
1576	$ annee_ref,
1577	$ day_ini)
1578
1579	CALL VTe(VTinca)
1580	CALL VTb(VTphysiq)
1581	#endif
1582	END IF
1583	c
1584	c
1585	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1586	! Nouvelle initialisation pour le rayonnement RRTM
1587	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1588
1589	call iniradia(klon,klev,paprs(1,1:klev+1))
1590
1591	ENDIF
1592	!
1593	! ************** Fin de IF ( debut ) *************
1594	!
1595	!
1596	! Incrementer le compteur de la physique
1597	!
1598	itap = itap + 1
1599	julien = MOD(NINT(xjour),360)
1600	if (julien .eq. 0) julien = 360
1601
1602	!
1603	! Update fraction of the sub-surfaces (pctsrf) and
1604	! initialize, where a new fraction has appeared, all variables depending
1605	! on the surface fraction.
1606	!
1607	CALL change_srf_frac(itap, dtime, julien,
1608	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1609
1610
1611	! Tendances bidons pour les processus qui n'affectent pas certaines
1612	! variables.
1613	du0(:,:)=0.
1614	dv0(:,:)=0.
1615	dq0(:,:)=0.
1616	dql0(:,:)=0.
1617	c
1618	c Mettre a zero des variables de sortie (pour securite)
1619	c
1620	DO i = 1, klon
1621	d_ps(i) = 0.0
1622	ENDDO
1623	DO k = 1, klev
1624	DO i = 1, klon
1625	d_t(i,k) = 0.0
1626	d_u(i,k) = 0.0
1627	d_v(i,k) = 0.0
1628	ENDDO
1629	ENDDO
1630	DO iq = 1, nqmax
1631	DO k = 1, klev
1632	DO i = 1, klon
1633	d_qx(i,k,iq) = 0.0
1634	ENDDO
1635	ENDDO
1636	ENDDO
1637	da(:,:)=0.
1638	mp(:,:)=0.
1639	phi(:,:,:)=0.
1640	c
1641	c Ne pas affecter les valeurs entrees de u, v, h, et q
1642	c
1643	DO k = 1, klev
1644	DO i = 1, klon
1645	t_seri(i,k) = t(i,k)
1646	u_seri(i,k) = u(i,k)
1647	v_seri(i,k) = v(i,k)
1648	q_seri(i,k) = qx(i,k,ivap)
1649	ql_seri(i,k) = qx(i,k,iliq)
1650	qs_seri(i,k) = 0.
1651	ENDDO
1652	ENDDO
1653	IF (nqmax.GE.3) THEN
1654	DO iq = 3, nqmax
1655	DO k = 1, klev
1656	DO i = 1, klon
1657	tr_seri(i,k,iq-2) = qx(i,k,iq)
1658	ENDDO
1659	ENDDO
1660	ENDDO
1661	ELSE
1662	DO k = 1, klev
1663	DO i = 1, klon
1664	tr_seri(i,k,1) = 0.0
1665	ENDDO
1666	ENDDO
1667	ENDIF
1668	C
1669	DO i = 1, klon
1670	ztsol(i) = 0.
1671	ENDDO
1672	DO nsrf = 1, nbsrf
1673	DO i = 1, klon
1674	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1675	ENDDO
1676	ENDDO
1677	cIM
1678	IF (ip_ebil_phy.ge.1) THEN
1679	ztit='after dynamic'
1680	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1681	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1682	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1683	C Comme les tendances de la physique sont ajoute dans la dynamique,
1684	C on devrait avoir que la variation d'entalpie par la dynamique
1685	C est egale a la variation de la physique au pas de temps precedent.
1686	C Donc la somme de ces 2 variations devrait etre nulle.
1687	call diagphy(airephy,ztit,ip_ebil_phy
1688	e , zero_v, zero_v, zero_v, zero_v, zero_v
1689	e , zero_v, zero_v, zero_v, ztsol
1690	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1691	s , fs_bound, fq_bound )
1692	END IF
1693
1694	c Diagnostiquer la tendance dynamique
1695	c
1696	IF (ancien_ok) THEN
1697	DO k = 1, klev
1698	DO i = 1, klon
1699	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1700	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1701	ENDDO
1702	ENDDO
1703	ELSE
1704	DO k = 1, klev
1705	DO i = 1, klon
1706	d_t_dyn(i,k) = 0.0
1707	d_q_dyn(i,k) = 0.0
1708	ENDDO
1709	ENDDO
1710	ancien_ok = .TRUE.
1711	ENDIF
1712	c
1713	c Ajouter le geopotentiel du sol:
1714	c
1715	DO k = 1, klev
1716	DO i = 1, klon
1717	zphi(i,k) = pphi(i,k) + pphis(i)
1718	ENDDO
1719	ENDDO
1720	c
1721	c Verifier les temperatures
1722	c
1723	cIM BEG
1724	IF (check) THEN
1725	amn=MIN(ftsol(1,is_ter),1000.)
1726	amx=MAX(ftsol(1,is_ter),-1000.)
1727	DO i=2, klon
1728	amn=MIN(ftsol(i,is_ter),amn)
1729	amx=MAX(ftsol(i,is_ter),amx)
1730	ENDDO
1731	c
1732	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1733	ENDIF !(check) THEN
1734	cIM END
1735	c
1736	CALL hgardfou(t_seri,ftsol,'debutphy')
1737	c
1738	cIM BEG
1739	IF (check) THEN
1740	amn=MIN(ftsol(1,is_ter),1000.)
1741	amx=MAX(ftsol(1,is_ter),-1000.)
1742	DO i=2, klon
1743	amn=MIN(ftsol(i,is_ter),amn)
1744	amx=MAX(ftsol(i,is_ter),amx)
1745	ENDDO
1746	c
1747	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1748	ENDIF !(check) THEN
1749	cIM END
1750	c
1751	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1752	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1753	c
1754	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1755	WRITE(lunout,*)' PHYS cond julien ',julien
1756	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1757	ENDIF
1758	c
1759	c Re-evaporer l'eau liquide nuageuse
1760	c
1761	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1762	DO i = 1, klon
1763	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1764	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1765	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1766	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1767	zb = MAX(0.0,ql_seri(i,k))
1768	za = - MAX(0.0,ql_seri(i,k))
1769	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1770	t_seri(i,k) = t_seri(i,k) + za
1771	q_seri(i,k) = q_seri(i,k) + zb
1772	ql_seri(i,k) = 0.0
1773	d_t_eva(i,k) = za
1774	d_q_eva(i,k) = zb
1775	ENDDO
1776	ENDDO
1777	cIM
1778	IF (ip_ebil_phy.ge.2) THEN
1779	ztit='after reevap'
1780	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1781	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1782	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1783	call diagphy(airephy,ztit,ip_ebil_phy
1784	e , zero_v, zero_v, zero_v, zero_v, zero_v
1785	e , zero_v, zero_v, zero_v, ztsol
1786	e , d_h_vcol, d_qt, d_ec
1787	s , fs_bound, fq_bound )
1788	C
1789	END IF
1790
1791	c
1792	c=========================================================================
1793	! Calculs de l'orbite.
1794	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1795	! doit donc etre placé avant radlwsw et pbl_surface
1796
1797	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1798	! solarlong0
1799
1800	if (solarlong0<-999.) then
1801	CALL orbite(FLOAT(julien),zlongi,dist)
1802	else
1803	zlongi=solarlong0 ! longitude solaire vraie
1804	dist=1. ! distance au soleil / moyenne
1805	endif
1806
1807	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
1808
1809	! Avec ou sans cycle diurne
1810	IF (cycle_diurne) THEN
1811	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1812	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1813	ELSE
1814	rmu0 = -999.999
1815	ENDIF
1816
1817	if (mydebug) then
1818	call writefield_phy('u_seri',u_seri,llm)
1819	call writefield_phy('v_seri',v_seri,llm)
1820	call writefield_phy('t_seri',t_seri,llm)
1821	call writefield_phy('q_seri',q_seri,llm)
1822	endif
1823
1824	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1825	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1826	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1827	c turbulent du couche limit.
1828	c
1829	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1830	c ecriture des fihiers hist_XXXX.nc, ces sont :
1831	c qsol, zq2m, s_pblh, s_lcl,
1832	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1833	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1834	c zxrugs, zu10m, zv10m, fder,
1835	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1836	c frugs, agesno, fsollw, fsolsw,
1837	c d_ts, fevap, fluxlat, t2m,
1838	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1839	c
1840	c Certains ne sont pas utiliser du tout :
1841	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1842	c
1843
1844	CALL pbl_surface(
1845	e dtime, date0, itap, julien,
1846	e debut, lafin,
1847	e rlon, rlat, rugoro, rmu0,
1848	e rain_fall, snow_fall, solsw, sollw,
1849	e t_seri, q_seri, u_seri, v_seri,
1850	e pplay, paprs, pctsrf,
1851	+ ftsol, falb1, falb2, u10m, v10m,
1852	s sollwdown, cdragh, cdragm, yu1, yv1,
1853	s albsol1, albsol2, sens, evap,
1854	s zxtsol, zxfluxlat, zt2m, qsat2m,
1855	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1856	s ycoefh, slab_wfbils,
1857	d qsol, zq2m, s_pblh, s_lcl,
1858	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1859	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1860	d zxrugs, zu10m, zv10m, fder,
1861	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1862	d frugs, agesno, fsollw, fsolsw,
1863	d d_ts, fevap, fluxlat, t2m,
1864	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1865	- dsens, devap, zxsnow,
1866	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1867
1868
1869	!-----------------------------------------------------------------------------------------
1870	! ajout des tendances de la diffusion turbulente
1871	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1872	!-----------------------------------------------------------------------------------------
1873
1874	if (mydebug) then
1875	call writefield_phy('u_seri',u_seri,llm)
1876	call writefield_phy('v_seri',v_seri,llm)
1877	call writefield_phy('t_seri',t_seri,llm)
1878	call writefield_phy('q_seri',q_seri,llm)
1879	endif
1880
1881
1882	IF (ip_ebil_phy.ge.2) THEN
1883	ztit='after surface_main'
1884	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1885	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1886	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1887	call diagphy(airephy,ztit,ip_ebil_phy
1888	e , zero_v, zero_v, zero_v, zero_v, sens
1889	e , evap , zero_v, zero_v, ztsol
1890	e , d_h_vcol, d_qt, d_ec
1891	s , fs_bound, fq_bound )
1892	END IF
1893
1894	c =================================================================== c
1895	c Calcul de Qsat
1896
1897	DO k = 1, klev
1898	DO i = 1, klon
1899	zx_t = t_seri(i,k)
1900	IF (thermcep) THEN
1901	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1902	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1903	zx_qs = MIN(0.5,zx_qs)
1904	zcor = 1./(1.-retv*zx_qs)
1905	zx_qs = zx_qs*zcor
1906	ELSE
1907	IF (zx_t.LT.t_coup) THEN
1908	zx_qs = qsats(zx_t)/pplay(i,k)
1909	ELSE
1910	zx_qs = qsatl(zx_t)/pplay(i,k)
1911	ENDIF
1912	ENDIF
1913	zqsat(i,k)=zx_qs
1914	ENDDO
1915	ENDDO
1916
1917	if (prt_level.ge.1) then
1918	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1919	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1920	endif
1921	c
1922	c Appeler la convection (au choix)
1923	c
1924	DO k = 1, klev
1925	DO i = 1, klon
1926	conv_q(i,k) = d_q_dyn(i,k)
1927	. + d_q_vdf(i,k)/dtime
1928	conv_t(i,k) = d_t_dyn(i,k)
1929	. + d_t_vdf(i,k)/dtime
1930	ENDDO
1931	ENDDO
1932	IF (check) THEN
1933	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1934	WRITE(lunout,*) "avantcon=", za
1935	ENDIF
1936	zx_ajustq = .FALSE.
1937	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1938	IF (zx_ajustq) THEN
1939	DO i = 1, klon
1940	z_avant(i) = 0.0
1941	ENDDO
1942	DO k = 1, klev
1943	DO i = 1, klon
1944	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1945	. *(paprs(i,k)-paprs(i,k+1))/RG
1946	ENDDO
1947	ENDDO
1948	ENDIF
1949
1950	c Calcule de vitesse verticale a partir de flux de masse verticale
1951	DO k = 1, klev
1952	DO i = 1, klon
1953	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1954	END DO
1955	END DO
1956
1957	IF (iflag_con.EQ.1) THEN
1958	stop'reactiver le call conlmd dans physiq.F'
1959	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1960	c . d_t_con, d_q_con,
1961	c . rain_con, snow_con, ibas_con, itop_con)
1962	ELSE IF (iflag_con.EQ.2) THEN
1963	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1964	e conv_t, conv_q, -evap, omega,
1965	s d_t_con, d_q_con, rain_con, snow_con,
1966	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1967	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1968	d_u_con = 0.
1969	d_v_con = 0.
1970
1971	WHERE (rain_con < 0.) rain_con = 0.
1972	WHERE (snow_con < 0.) snow_con = 0.
1973	DO i = 1, klon
1974	ibas_con(i) = klev+1 - kcbot(i)
1975	itop_con(i) = klev+1 - kctop(i)
1976	ENDDO
1977	ELSE IF (iflag_con.GE.3) THEN
1978	c nb of tracers for the KE convection:
1979	c MAF la partie traceurs est faite dans phytrac
1980	c on met ntra=1 pour limiter les appels mais on peut
1981	c supprimer les calculs / ftra.
1982	ntra = 1
1983
1984	c=====================================================================================
1985	cajout pour la parametrisation des poches froides:
1986	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1987	do k=1,klev
1988	do i=1,klon
1989	if (iflag_wake.eq.1) then
1990	t_wake(i,k) = t_seri(i,k)
1991	. +(1-wake_s(i))*wake_deltat(i,k)
1992	q_wake(i,k) = q_seri(i,k)
1993	. +(1-wake_s(i))*wake_deltaq(i,k)
1994	t_undi(i,k) = t_seri(i,k)
1995	. -wake_s(i)*wake_deltat(i,k)
1996	q_undi(i,k) = q_seri(i,k)
1997	. -wake_s(i)*wake_deltaq(i,k)
1998	else
1999	t_wake(i,k) = t_seri(i,k)
2000	q_wake(i,k) = q_seri(i,k)
2001	t_undi(i,k) = t_seri(i,k)
2002	q_undi(i,k) = q_seri(i,k)
2003	endif
2004	enddo
2005	enddo
2006
2007	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2008	cc-- pour le soulevement des particules dans le modele convectif
2009	c
2010	do i = 1,klon
2011	ALE(i) = 0.
2012	ALP(i) = 0.
2013	enddo
2014	c
2015	ccalcul de ale_wake et alp_wake
2016	do i = 1,klon
2017	if (iflag_wake.eq.1) then
2018	ale_wake(i) = 0.5wake_cstar(i)*2
2019	alp_wake(i) = wake_fip(i)
2020	else
2021	ale_wake(i) = 0.
2022	alp_wake(i) = 0.
2023	endif
2024	enddo
2025	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2026	cdans le thermique sinon
2027	if (iflag_coupl.eq.0) then
2028	if (debut) print*,'ALE et ALP imposes'
2029	do i = 1,klon
2030	con ne couple que ale
2031	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2032	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2033	con ne couple que alp
2034	c ALP(i) = alp_wake(i) + Alp_bl(i)
2035	ALP(i) = alp_wake(i) + alp_bl_prescr
2036	enddo
2037	else
2038	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2039	do i = 1,klon
2040	ALE(i) = max(ale_wake(i),Ale_bl(i))
2041	ALP(i) = alp_wake(i) + Alp_bl(i)
2042	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2043	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2044	enddo
2045	endif
2046	do i=1,klon
2047	if (alp(i)>alp_max) then
2048	print*,'WARNING SUPER ALP (seuil=',alp_max,
2049	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2050	alp(i)=alp_max
2051	endif
2052	if (ale(i)>ale_max) then
2053	print*,'WARNING SUPER ALE (seuil=',ale_max,
2054	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2055	ale(i)=ale_max
2056	endif
2057	enddo
2058
2059	cfin calcul ale et alp
2060	c=================================================================================================
2061
2062
2063	c sb, oct02:
2064	c Schema de convection modularise et vectorise:
2065	c (driver commun aux versions 3 et 4)
2066	c
2067	IF (ok_cvl) THEN ! new driver for convectL
2068
2069	CALL concvl (iflag_con,iflag_clos,
2070	. dtime,paprs,pplay,t_undi,q_undi,
2071	. t_wake,q_wake,
2072	. u_seri,v_seri,tr_seri,nbtr,
2073	. ALE,ALP,
2074	. ema_work1,ema_work2,
2075	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2076	. rain_con, snow_con, ibas_con, itop_con, sigd,
2077	. upwd,dnwd,dnwd0,
2078	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2079	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2080	. pmflxr,pmflxs,da,phi,mp,
2081	. ftd,fqd,lalim_conv,wght_th)
2082
2083	cIM begin
2084	print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2085	.dnwd0(1,1),ftd(1,1),fqd(1,1)
2086	cIM end
2087	cIM cf. FH
2088	clwcon0=qcondc
2089	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2090
2091	ELSE ! ok_cvl
2092	c MAF conema3 ne contient pas les traceurs
2093	CALL conema3 (dtime,
2094	. paprs,pplay,t_seri,q_seri,
2095	. u_seri,v_seri,tr_seri,ntra,
2096	. ema_work1,ema_work2,
2097	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2098	. rain_con, snow_con, ibas_con, itop_con,
2099	. upwd,dnwd,dnwd0,bas,top,
2100	. Ma,cape,tvp,rflag,
2101	. pbase
2102	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2103	. ,clwcon0)
2104
2105	ENDIF ! ok_cvl
2106
2107	c
2108	c Correction precip
2109	rain_con = rain_con * cvl_corr
2110	snow_con = snow_con * cvl_corr
2111	c
2112
2113	IF (.NOT. ok_gust) THEN
2114	do i = 1, klon
2115	wd(i)=0.0
2116	enddo
2117	ENDIF
2118
2119	c =================================================================== c
2120	c Calcul des proprietes des nuages convectifs
2121	c
2122
2123	c calcul des proprietes des nuages convectifs
2124	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2125	call clouds_gno
2126	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2127
2128	c =================================================================== c
2129
2130	DO i = 1, klon
2131	ema_pcb(i) = pbase(i)
2132	ENDDO
2133	DO i = 1, klon
2134
2135	! L'idicage de itop_con peut cacher un pb potentiel
2136	! FH sous la dictee de JYG, CR
2137	ema_pct(i) = paprs(i,itop_con(i)+1)
2138
2139	if (itop_con(i).gt.klev-3) then
2140	print*,'La convection monte trop haut '
2141	print*,'itop_con(,',i,',)=',itop_con(i)
2142	endif
2143	ENDDO
2144	DO i = 1, klon
2145	ema_cbmf(i) = ema_workcbmf(i)
2146	ENDDO
2147	ELSE IF (iflag_con.eq.0) THEN
2148	write(lunout,*) 'On n appelle pas la convection'
2149	clwcon0=0.
2150	rnebcon0=0.
2151	d_t_con=0.
2152	d_q_con=0.
2153	d_u_con=0.
2154	d_v_con=0.
2155	rain_con=0.
2156	snow_con=0.
2157	bas=1
2158	top=1
2159	ELSE
2160	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2161	CALL abort
2162	ENDIF
2163
2164	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2165	c . d_u_con, d_v_con)
2166
2167	!-----------------------------------------------------------------------------------------
2168	! ajout des tendances de la diffusion turbulente
2169	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2170	!-----------------------------------------------------------------------------------------
2171
2172	if (mydebug) then
2173	call writefield_phy('u_seri',u_seri,llm)
2174	call writefield_phy('v_seri',v_seri,llm)
2175	call writefield_phy('t_seri',t_seri,llm)
2176	call writefield_phy('q_seri',q_seri,llm)
2177	endif
2178
2179	cIM
2180	IF (ip_ebil_phy.ge.2) THEN
2181	ztit='after convect'
2182	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2183	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2184	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2185	call diagphy(airephy,ztit,ip_ebil_phy
2186	e , zero_v, zero_v, zero_v, zero_v, zero_v
2187	e , zero_v, rain_con, snow_con, ztsol
2188	e , d_h_vcol, d_qt, d_ec
2189	s , fs_bound, fq_bound )
2190	END IF
2191	C
2192	IF (check) THEN
2193	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2194	WRITE(lunout,*)"aprescon=", za
2195	zx_t = 0.0
2196	za = 0.0
2197	DO i = 1, klon
2198	za = za + airephy(i)/FLOAT(klon)
2199	zx_t = zx_t + (rain_con(i)+
2200	. snow_con(i))*airephy(i)/FLOAT(klon)
2201	ENDDO
2202	zx_t = zx_t/za*dtime
2203	WRITE(lunout,*)"Precip=", zx_t
2204	ENDIF
2205	IF (zx_ajustq) THEN
2206	DO i = 1, klon
2207	z_apres(i) = 0.0
2208	ENDDO
2209	DO k = 1, klev
2210	DO i = 1, klon
2211	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2212	. *(paprs(i,k)-paprs(i,k+1))/RG
2213	ENDDO
2214	ENDDO
2215	DO i = 1, klon
2216	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2217	. /z_apres(i)
2218	ENDDO
2219	DO k = 1, klev
2220	DO i = 1, klon
2221	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2222	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2223	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2224	ENDIF
2225	ENDDO
2226	ENDDO
2227	ENDIF
2228	zx_ajustq=.FALSE.
2229
2230	c
2231	c=============================================================================
2232	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2233	cpour la couche limite diffuse pour l instant
2234	c
2235	if (iflag_wake.eq.1) then
2236	DO k=1,klev
2237	DO i=1,klon
2238	dt_dwn(i,k) = ftd(i,k)
2239	wdt_PBL(i,k) = 0.
2240	dq_dwn(i,k) = fqd(i,k)
2241	wdq_PBL(i,k) = 0.
2242	M_dwn(i,k) = dnwd0(i,k)
2243	M_up(i,k) = upwd(i,k)
2244	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2245	udt_PBL(i,k) = 0.
2246	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2247	udq_PBL(i,k) = 0.
2248	ENDDO
2249	ENDDO
2250	c
2251	ccalcul caracteristiques de la poche froide
2252	call calWAKE (paprs,pplay,dtime
2253	: ,t_seri,q_seri,omega
2254	: ,dt_dwn,dq_dwn,M_dwn,M_up
2255	: ,dt_a,dq_a,sigd
2256	: ,wdt_PBL,wdq_PBL
2257	: ,udt_PBL,udq_PBL
2258	o ,wake_deltat,wake_deltaq,wake_dth
2259	o ,wake_h,wake_s,wake_dens
2260	o ,wake_pe,wake_fip,wake_gfl
2261	o ,dt_wake,dq_wake
2262	o ,wake_k, t_undi,q_undi
2263	o ,wake_omgbdth,wake_dp_omgb
2264	o ,wake_dtKE,wake_dqKE
2265	o ,wake_dtPBL,wake_dqPBL
2266	o ,wake_omg,wake_dp_deltomg
2267	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2268	o ,wake_ddeltat,wake_ddeltaq)
2269	c
2270	!-----------------------------------------------------------------------------------------
2271	! ajout des tendances des poches froides
2272	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2273	! coherent avec les autres d_t_...
2274	d_t_wake(:,:)=dt_wake(:,:)*dtime
2275	d_q_wake(:,:)=dq_wake(:,:)*dtime
2276	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2277	!-----------------------------------------------------------------------------------------
2278
2279	endif
2280	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2281	c
2282	c===================================================================
2283	c Convection seche (thermiques ou ajustement)
2284	c===================================================================
2285	c
2286	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2287	s ,seuil_inversion,weak_inversion,dthmin)
2288
2289
2290
2291	d_t_ajsb(:,:)=0.
2292	d_q_ajsb(:,:)=0.
2293	d_t_ajs(:,:)=0.
2294	d_u_ajs(:,:)=0.
2295	d_v_ajs(:,:)=0.
2296	d_q_ajs(:,:)=0.
2297	clwcon0th(:,:)=0.
2298	c
2299	fm_therm(:,:)=0.
2300	entr_therm(:,:)=0.
2301	detr_therm(:,:)=0.
2302	c
2303	IF(prt_level>9)WRITE(lunout,*)
2304	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2305	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2306	if(iflag_thermals.lt.0) then
2307	c Rien
2308	c ====
2309	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2310
2311
2312	else
2313
2314	c Thermiques
2315	c ==========
2316	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2317	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2318
2319
2320	if (iflag_thermals.gt.1) then
2321	call calltherm(pdtphys
2322	s ,pplay,paprs,pphi,weak_inversion
2323	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2324	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2325	s ,fm_therm,entr_therm,detr_therm
2326	s ,zqasc,clwcon0th,lmax_th,ratqscth
2327	s ,ratqsdiff,zqsatth
2328	con rajoute ale et alp, et les caracteristiques de la couche alim
2329	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2330	endif
2331
2332
2333	c Ajustement sec
2334	c ==============
2335
2336	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2337	! a partir du sommet des thermiques.
2338	! Dans le cas contraire, on demarre au niveau 1.
2339
2340	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2341
2342	if(iflag_thermals.eq.0) then
2343	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2344	limbas(:)=1
2345	else
2346	limbas(:)=lmax_th(:)
2347	endif
2348
2349	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2350	! pour des test de convergence numerique.
2351	! Le nouveau ajsec est a priori mieux, meme pour le cas
2352	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2353	! non nulles numeriquement pour des mailles non concernees.
2354
2355	if (iflag_thermals.eq.0) then
2356	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2357	s , d_t_ajsb, d_q_ajsb)
2358	else
2359	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2360	s , d_t_ajsb, d_q_ajsb)
2361	endif
2362
2363	!-----------------------------------------------------------------------------------------
2364	! ajout des tendances de l'ajustement sec ou des thermiques
2365	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2366	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2367	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2368
2369	!-----------------------------------------------------------------------------------------
2370
2371	endif
2372
2373	endif
2374	c
2375	c===================================================================
2376	cIM
2377	IF (ip_ebil_phy.ge.2) THEN
2378	ztit='after dry_adjust'
2379	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2380	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2381	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2382	END IF
2383
2384
2385	c-------------------------------------------------------------------------
2386	c Caclul des ratqs
2387	c-------------------------------------------------------------------------
2388
2389	c print*,'calcul des ratqs'
2390	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2391	c ----------------
2392	c on ecrase le tableau ratqsc calcule par clouds_gno
2393	if (iflag_cldcon.eq.1) then
2394	do k=1,klev
2395	do i=1,klon
2396	if(ptconv(i,k)) then
2397	ratqsc(i,k)=ratqsbas
2398	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2399	else
2400	ratqsc(i,k)=0.
2401	endif
2402	enddo
2403	enddo
2404
2405	c-----------------------------------------------------------------------
2406	c par nversion de la fonction log normale
2407	c-----------------------------------------------------------------------
2408	else if (iflag_cldcon.eq.4) then
2409	ptconvth(:,:)=.false.
2410	ratqsc(:,:)=0.
2411	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2412	call clouds_gno
2413	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2414	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2415
2416	c-----------------------------------------------------------------------
2417	c par calcul direct de l'ecart-type
2418	c-----------------------------------------------------------------------
2419
2420	else if (iflag_cldcon>=5) then
2421	wmax_th(:)=0.
2422	zmax_th(:)=0.
2423	do k=1,klev
2424	do i=1,klon
2425	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2426	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2427	enddo
2428	enddo
2429	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2430	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
2431
2432	c On impose que l'air autour de la fraction couverte par le thermique
2433	c plus son air detraine durant tau_overturning_th soit superieur
2434	c a 0.1 q_seri
2435	zz=0.1
2436	do k=1,klev
2437	do i=1,klon
2438	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2439	s tau_overturning_th(i)*detr_therm(i,k)
2440	s *rg/(paprs(i,k)-paprs(i,k+1))
2441	znum=(1.-zz)*q_seri(i,k)
2442	zden=zqasc(i,k)-zz*q_seri(i,k)
2443	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2444	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2445	enddo
2446	enddo
2447
2448	if(iflag_cldcon==5) then
2449	do k=1,klev
2450	do i=1,klon
2451	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2452	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2453	enddo
2454	enddo
2455	else if(iflag_cldcon==6) then
2456	do k=1,klev
2457	do i=1,klon
2458	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2459	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2460	enddo
2461	enddo
2462	endif
2463
2464	endif
2465
2466	c ratqs stables
2467	c -------------
2468
2469	if (iflag_ratqs.eq.0) then
2470
2471	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2472	do k=1,klev
2473	do i=1, klon
2474	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2475	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2476	enddo
2477	enddo
2478
2479	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2480	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2481	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2482	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2483	! Il s'agit de differents tests dans la phase de reglage du modele
2484	! avec thermiques.
2485
2486	else if (iflag_ratqs.eq.1) then
2487
2488	do k=1,klev
2489	do i=1, klon
2490	if (pplay(i,k).ge.60000.) then
2491	ratqss(i,k)=ratqsbas
2492	else if ((pplay(i,k).ge.30000.).and.
2493	s (pplay(i,k).lt.60000.)) then
2494	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2495	s (60000.-pplay(i,k))/(60000.-30000.)
2496	else
2497	ratqss(i,k)=ratqshaut
2498	endif
2499	enddo
2500	enddo
2501
2502	else
2503
2504	do k=1,klev
2505	do i=1, klon
2506	if (pplay(i,k).ge.60000.) then
2507	ratqss(i,k)=ratqsbas
2508	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2509	else if ((pplay(i,k).ge.30000.).and.
2510	s (pplay(i,k).lt.60000.)) then
2511	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2512	s (60000.-pplay(i,k))/(60000.-30000.)
2513	else
2514	ratqss(i,k)=ratqshaut
2515	endif
2516	enddo
2517	enddo
2518	endif
2519
2520
2521
2522
2523	c ratqs final
2524	c -----------
2525
2526	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2527	s .or.iflag_cldcon.ge.4) then
2528
2529	! On ajoute une constante au ratqsc*2 pour tenir compte de
2530	! fluctuations turbulentes de petite echelle
2531
2532	do k=1,klev
2533	do i=1,klon
2534	if ((fm_therm(i,k).gt.1.e-10)) then
2535	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2536	endif
2537	enddo
2538	enddo
2539
2540	! les ratqs sont une conbinaison de ratqss et ratqsc
2541	! 1800s, en dur pour le moment, est le temps de
2542	! relaxation des ratqs
2543
2544	facteur=exp(-pdtphys/1800.)
2545
2546	print*,'WARNING ratqs a revoir '
2547	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2548	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2549
2550	else
2551	! on ne prend que le ratqs stable pour fisrtilp
2552	ratqs(:,:)=ratqss(:,:)
2553	endif
2554
2555
2556	c
2557	c Appeler le processus de condensation a grande echelle
2558	c et le processus de precipitation
2559	c-------------------------------------------------------------------------
2560	CALL fisrtilp(dtime,paprs,pplay,
2561	. t_seri, q_seri,ptconv,ratqs,
2562	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2563	. rain_lsc, snow_lsc,
2564	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2565	. frac_impa, frac_nucl,
2566	. prfl, psfl, rhcl)
2567
2568	WHERE (rain_lsc < 0) rain_lsc = 0.
2569	WHERE (snow_lsc < 0) snow_lsc = 0.
2570	!-----------------------------------------------------------------------------------------
2571	! ajout des tendances de la diffusion turbulente
2572	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2573	!-----------------------------------------------------------------------------------------
2574	DO k = 1, klev
2575	DO i = 1, klon
2576	cldfra(i,k) = rneb(i,k)
2577	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2578	ENDDO
2579	ENDDO
2580	IF (check) THEN
2581	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2582	WRITE(lunout,*)"apresilp=", za
2583	zx_t = 0.0
2584	za = 0.0
2585	DO i = 1, klon
2586	za = za + airephy(i)/FLOAT(klon)
2587	zx_t = zx_t + (rain_lsc(i)
2588	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2589	ENDDO
2590	zx_t = zx_t/za*dtime
2591	WRITE(lunout,*)"Precip=", zx_t
2592	ENDIF
2593	cIM
2594	IF (ip_ebil_phy.ge.2) THEN
2595	ztit='after fisrt'
2596	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2597	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2598	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2599	call diagphy(airephy,ztit,ip_ebil_phy
2600	e , zero_v, zero_v, zero_v, zero_v, zero_v
2601	e , zero_v, rain_lsc, snow_lsc, ztsol
2602	e , d_h_vcol, d_qt, d_ec
2603	s , fs_bound, fq_bound )
2604	END IF
2605
2606	if (mydebug) then
2607	call writefield_phy('u_seri',u_seri,llm)
2608	call writefield_phy('v_seri',v_seri,llm)
2609	call writefield_phy('t_seri',t_seri,llm)
2610	call writefield_phy('q_seri',q_seri,llm)
2611	endif
2612
2613	c
2614	c-------------------------------------------------------------------
2615	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2616	c-------------------------------------------------------------------
2617
2618	c 1. NUAGES CONVECTIFS
2619	c
2620	cIM cf FH
2621	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2622	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2623	snow_tiedtke=0.
2624	c print*,'avant calcul de la pseudo precip '
2625	c print*,'iflag_cldcon',iflag_cldcon
2626	if (iflag_cldcon.eq.-1) then
2627	rain_tiedtke=rain_con
2628	else
2629	c print*,'calcul de la pseudo precip '
2630	rain_tiedtke=0.
2631	c print*,'calcul de la pseudo precip 0'
2632	do k=1,klev
2633	do i=1,klon
2634	if (d_q_con(i,k).lt.0.) then
2635	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2636	s *(paprs(i,k)-paprs(i,k+1))/rg
2637	endif
2638	enddo
2639	enddo
2640	endif
2641	c
2642	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2643	c
2644
2645	c Nuages diagnostiques pour Tiedtke
2646	CALL diagcld1(paprs,pplay,
2647	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2648	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2649	. diafra,dialiq)
2650	DO k = 1, klev
2651	DO i = 1, klon
2652	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2653	cldliq(i,k) = dialiq(i,k)
2654	cldfra(i,k) = diafra(i,k)
2655	ENDIF
2656	ENDDO
2657	ENDDO
2658
2659	ELSE IF (iflag_cldcon.ge.3) THEN
2660	c On prend pour les nuages convectifs le max du calcul de la
2661	c convection et du calcul du pas de temps precedent diminue d'un facteur
2662	c facttemps
2663	facteur = pdtphys *facttemps
2664	do k=1,klev
2665	do i=1,klon
2666	rnebcon(i,k)=rnebcon(i,k)*facteur
2667	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2668	s then
2669	rnebcon(i,k)=rnebcon0(i,k)
2670	clwcon(i,k)=clwcon0(i,k)
2671	endif
2672	enddo
2673	enddo
2674
2675	c
2676	cjq - introduce the aerosol direct and first indirect radiative forcings
2677	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2678	IF (ok_ade.OR.ok_aie) THEN
2679	IF ( .NOT. aerosol_couple ) THEN
2680	! Get sulfate aerosol distribution
2681	CALL readsulfate(rjourvrai, debut, sulfate)
2682	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2683
2684	! Calculate aerosol optical properties (Olivier Boucher)
2685	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2686	. tau_ae, piz_ae, cg_ae, aerindex)
2687	ENDIF
2688	ELSE
2689	tau_ae(:,:,:)=0.0
2690	piz_ae(:,:,:)=0.0
2691	cg_ae(:,:,:)=0.0
2692	ENDIF
2693
2694	cIM calcul nuages par le simulateur ISCCP
2695	c
2696	#ifdef histISCCP
2697	IF (ok_isccp) THEN
2698	c
2699	c
2700	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2701	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2702	#include "calcul_simulISCCP.h"
2703	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2704	ENDIF !ok_isccp
2705	#endif
2706
2707	c On prend la somme des fractions nuageuses et des contenus en eau
2708	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2709	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2710
2711	ENDIF
2712	c
2713	c 2. NUAGES STARTIFORMES
2714	c
2715	IF (ok_stratus) THEN
2716	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2717	DO k = 1, klev
2718	DO i = 1, klon
2719	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2720	cldliq(i,k) = dialiq(i,k)
2721	cldfra(i,k) = diafra(i,k)
2722	ENDIF
2723	ENDDO
2724	ENDDO
2725	ENDIF
2726	c
2727	c Precipitation totale
2728	c
2729	DO i = 1, klon
2730	rain_fall(i) = rain_con(i) + rain_lsc(i)
2731	snow_fall(i) = snow_con(i) + snow_lsc(i)
2732	ENDDO
2733	cIM
2734	IF (ip_ebil_phy.ge.2) THEN
2735	ztit="after diagcld"
2736	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2737	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2738	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2739	END IF
2740	c
2741	c Calculer l'humidite relative pour diagnostique
2742	c
2743	DO k = 1, klev
2744	DO i = 1, klon
2745	zx_t = t_seri(i,k)
2746	IF (thermcep) THEN
2747	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2748	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2749	zx_qs = MIN(0.5,zx_qs)
2750	zcor = 1./(1.-retv*zx_qs)
2751	zx_qs = zx_qs*zcor
2752	ELSE
2753	IF (zx_t.LT.t_coup) THEN
2754	zx_qs = qsats(zx_t)/pplay(i,k)
2755	ELSE
2756	zx_qs = qsatl(zx_t)/pplay(i,k)
2757	ENDIF
2758	ENDIF
2759	zx_rh(i,k) = q_seri(i,k)/zx_qs
2760	zqsat(i,k)=zx_qs
2761	ENDDO
2762	ENDDO
2763
2764	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2765	c equivalente a 2m (tpote) pour diagnostique
2766	c
2767	DO i = 1, klon
2768	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2769	IF (thermcep) THEN
2770	IF(zt2m(i).LT.RTT) then
2771	Lheat=RLSTT
2772	ELSE
2773	Lheat=RLVTT
2774	ENDIF
2775	ELSE
2776	IF (zt2m(i).LT.RTT) THEN
2777	Lheat=RLSTT
2778	ELSE
2779	Lheat=RLVTT
2780	ENDIF
2781	ENDIF
2782	tpote(i) = tpot(i)*
2783	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2784	ENDDO
2785
2786	IF (config_inca /= 'none') THEN
2787	#ifdef INCA
2788	CALL VTe(VTphysiq)
2789	CALL VTb(VTinca)
2790	calday = FLOAT(julien) + gmtime
2791
2792	IF (config_inca == 'aero') THEN
2793	CALL AEROSOL_METEO_CALC(
2794	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2795	$ prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2796	END IF
2797
2798	zxsnow_dummy(:) = 0.0
2799
2800	CALL chemhook_begin (calday,
2801	$ julien,
2802	$ gmtime,
2803	$ pctsrf(1,1),
2804	$ rlat,
2805	$ rlon,
2806	$ airephy,
2807	$ paprs,
2808	$ pplay,
2809	$ ycoefh,
2810	$ pphi,
2811	$ t_seri,
2812	$ u,
2813	$ v,
2814	$ wo,
2815	$ q_seri,
2816	$ zxtsol,
2817	$ zxsnow_dummy,
2818	$ solsw,
2819	$ albsol1,
2820	$ rain_fall,
2821	$ snow_fall,
2822	$ itop_con,
2823	$ ibas_con,
2824	$ cldfra,
2825	$ iim,
2826	$ jjm,
2827	$ tr_seri,
2828	$ ftsol,
2829	$ paprs,
2830	$ cdragh,
2831	$ cdragm,
2832	$ pctsrf,
2833	$ pdtphys,
2834	$ itap)
2835
2836	CALL VTe(VTinca)
2837	CALL VTb(VTphysiq)
2838	#endif
2839	END IF !config_inca /= 'none'
2840	c
2841	c Calculer les parametres optiques des nuages et quelques
2842	c parametres pour diagnostiques:
2843	c
2844
2845	IF (aerosol_couple) THEN
2846	sulfate(:,:) = ccm(:,:,1)
2847	sulfate_pi(:,:) = ccm(:,:,2)
2848	ENDIF
2849
2850	if (ok_newmicro) then
2851	CALL newmicro (paprs, pplay,ok_newmicro,
2852	. t_seri, cldliq, cldfra, cldtau, cldemi,
2853	. cldh, cldl, cldm, cldt, cldq,
2854	. flwp, fiwp, flwc, fiwc,
2855	e ok_aie,
2856	e sulfate, sulfate_pi,
2857	e bl95_b0, bl95_b1,
2858	s cldtaupi, re, fl)
2859	else
2860	CALL nuage (paprs, pplay,
2861	. t_seri, cldliq, cldfra, cldtau, cldemi,
2862	. cldh, cldl, cldm, cldt, cldq,
2863	e ok_aie,
2864	e sulfate, sulfate_pi,
2865	e bl95_b0, bl95_b1,
2866	s cldtaupi, re, fl)
2867
2868	endif
2869	c
2870	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2871	c
2872	IF (MOD(itaprad,radpas).EQ.0) THEN
2873
2874	DO i = 1, klon
2875	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2876	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2877	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2878	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2879	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2880	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2881	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2882	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2883	ENDDO
2884
2885	if (mydebug) then
2886	call writefield_phy('u_seri',u_seri,llm)
2887	call writefield_phy('v_seri',v_seri,llm)
2888	call writefield_phy('t_seri',t_seri,llm)
2889	call writefield_phy('q_seri',q_seri,llm)
2890	endif
2891
2892	IF (aerosol_couple) THEN
2893	#ifdef INCA
2894	CALL radlwsw_inca
2895	e (kdlon,kflev,dist, rmu0, fract, solaire,
2896	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2897	e wo,
2898	e cldfra, cldemi, cldtau,
2899	s heat,heat0,cool,cool0,radsol,albpla,
2900	s topsw,toplw,solsw,sollw,
2901	s sollwdown,
2902	s topsw0,toplw0,solsw0,sollw0,
2903	s lwdn0, lwdn, lwup0, lwup,
2904	s swdn0, swdn, swup0, swup,
2905	e ok_ade, ok_aie,
2906	e tau_inca, piz_inca, cg_inca,
2907	s topswad_inca, solswad_inca,
2908	s topswad0_inca, solswad0_inca,
2909	s topsw_inca, topsw0_inca,
2910	s solsw_inca, solsw0_inca,
2911	e cldtaupi,
2912	s topswai_inca, solswai_inca)
2913	#endif
2914	ELSE
2915	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2916	e (dist, rmu0, fract,
2917	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2918	e wo,
2919	e cldfra, cldemi, cldtau,
2920	s heat,heat0,cool,cool0,radsol,albpla,
2921	s topsw,toplw,solsw,sollw,
2922	s sollwdown,
2923	s topsw0,toplw0,solsw0,sollw0,
2924	s lwdn0, lwdn, lwup0, lwup,
2925	s swdn0, swdn, swup0, swup,
2926	e ok_ade, ok_aie, ! new for aerosol radiative effects
2927	e tau_ae, piz_ae, cg_ae, ! ="=
2928	s topswad, solswad, ! ="=
2929	e cldtaupi, ! ="=
2930	s topswai, solswai,zqsat,flwc,fiwc) ! ="=
2931	ENDIF
2932	itaprad = 0
2933	ENDIF
2934	itaprad = itaprad + 1
2935
2936	if (iflag_radia.eq.0) then
2937	print *,'--------------------------------------------------'
2938	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2939	print *,'>>>> heat et cool mis a zero '
2940	print *,'--------------------------------------------------'
2941	heat=0.
2942	cool=0.
2943	endif
2944
2945
2946	c
2947	c Ajouter la tendance des rayonnements (tous les pas)
2948	c
2949	DO k = 1, klev
2950	DO i = 1, klon
2951	t_seri(i,k) = t_seri(i,k)
2952	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
2953	ENDDO
2954	ENDDO
2955	c
2956	if (mydebug) then
2957	call writefield_phy('u_seri',u_seri,llm)
2958	call writefield_phy('v_seri',v_seri,llm)
2959	call writefield_phy('t_seri',t_seri,llm)
2960	call writefield_phy('q_seri',q_seri,llm)
2961	endif
2962
2963	cIM
2964	IF (ip_ebil_phy.ge.2) THEN
2965	ztit='after rad'
2966	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2967	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2968	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2969	call diagphy(airephy,ztit,ip_ebil_phy
2970	e , topsw, toplw, solsw, sollw, zero_v
2971	e , zero_v, zero_v, zero_v, ztsol
2972	e , d_h_vcol, d_qt, d_ec
2973	s , fs_bound, fq_bound )
2974	END IF
2975	c
2976	c
2977	c Calculer l'hydrologie de la surface
2978	c
2979	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2980	c . agesno, ftsol,fqsurf,fsnow, ruis)
2981	c
2982
2983	c
2984	c Calculer le bilan du sol et la derive de temperature (couplage)
2985	c
2986	DO i = 1, klon
2987	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2988	c a la demande de JLD
2989	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2990	ENDDO
2991	c
2992	cmoddeblott(jan95)
2993	c Appeler le programme de parametrisation de l'orographie
2994	c a l'echelle sous-maille:
2995	c
2996	IF (ok_orodr) THEN
2997	c
2998	c selection des points pour lesquels le shema est actif:
2999	igwd=0
3000	DO i=1,klon
3001	itest(i)=0
3002	c IF ((zstd(i).gt.10.0)) THEN
3003	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3004	itest(i)=1
3005	igwd=igwd+1
3006	idx(igwd)=i
3007	ENDIF
3008	ENDDO
3009	c igwdim=MAX(1,igwd)
3010	c
3011	IF (ok_strato) THEN
3012
3013	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3014	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3015	e igwd,idx,itest,
3016	e t_seri, u_seri, v_seri,
3017	s zulow, zvlow, zustrdr, zvstrdr,
3018	s d_t_oro, d_u_oro, d_v_oro)
3019
3020	ELSE
3021	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3022	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3023	e igwd,idx,itest,
3024	e t_seri, u_seri, v_seri,
3025	s zulow, zvlow, zustrdr, zvstrdr,
3026	s d_t_oro, d_u_oro, d_v_oro)
3027	ENDIF
3028	c
3029	c ajout des tendances
3030	!-----------------------------------------------------------------------------------------
3031	! ajout des tendances de la trainee de l'orographie
3032	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3033	!-----------------------------------------------------------------------------------------
3034	c
3035	ENDIF ! fin de test sur ok_orodr
3036	c
3037	if (mydebug) then
3038	call writefield_phy('u_seri',u_seri,llm)
3039	call writefield_phy('v_seri',v_seri,llm)
3040	call writefield_phy('t_seri',t_seri,llm)
3041	call writefield_phy('q_seri',q_seri,llm)
3042	endif
3043
3044	IF (ok_orolf) THEN
3045	c
3046	c selection des points pour lesquels le shema est actif:
3047	igwd=0
3048	DO i=1,klon
3049	itest(i)=0
3050	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3051	itest(i)=1
3052	igwd=igwd+1
3053	idx(igwd)=i
3054	ENDIF
3055	ENDDO
3056	c igwdim=MAX(1,igwd)
3057	c
3058	IF (ok_strato) THEN
3059
3060	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3061	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3062	e igwd,idx,itest,
3063	e t_seri, u_seri, v_seri,
3064	s zulow, zvlow, zustrli, zvstrli,
3065	s d_t_lif, d_u_lif, d_v_lif )
3066
3067	ELSE
3068	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3069	e rlat,zmea,zstd,zpic,
3070	e itest,
3071	e t_seri, u_seri, v_seri,
3072	s zulow, zvlow, zustrli, zvstrli,
3073	s d_t_lif, d_u_lif, d_v_lif)
3074	ENDIF
3075	c
3076	!-----------------------------------------------------------------------------------------
3077	! ajout des tendances de la portance de l'orographie
3078	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3079	!-----------------------------------------------------------------------------------------
3080	c
3081	ENDIF ! fin de test sur ok_orolf
3082	C HINES GWD PARAMETRIZATION
3083
3084	IF (ok_hines) then
3085
3086	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3087	i rlat,t_seri,u_seri,v_seri,
3088	o zustrhi,zvstrhi,
3089	o d_t_hin, d_u_hin, d_v_hin)
3090	c
3091	c ajout des tendances
3092	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3093
3094	ENDIF
3095	c
3096
3097	c
3098	cIM cf. FLott BEG
3099	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3100
3101	if (mydebug) then
3102	call writefield_phy('u_seri',u_seri,llm)
3103	call writefield_phy('v_seri',v_seri,llm)
3104	call writefield_phy('t_seri',t_seri,llm)
3105	call writefield_phy('q_seri',q_seri,llm)
3106	endif
3107
3108	DO i = 1, klon
3109	zustrph(i)=0.
3110	zvstrph(i)=0.
3111	ENDDO
3112	DO k = 1, klev
3113	DO i = 1, klon
3114	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3115	c (paprs(i,k)-paprs(i,k+1))/rg
3116	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3117	c (paprs(i,k)-paprs(i,k+1))/rg
3118	ENDDO
3119	ENDDO
3120	c
3121	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3122	c
3123	IF (is_sequential) THEN
3124
3125	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3126	C ra,rg,romega,
3127	C rlat,rlon,pphis,
3128	C zustrdr,zustrli,zustrph,
3129	C zvstrdr,zvstrli,zvstrph,
3130	C paprs,u,v,
3131	C aam, torsfc)
3132	ENDIF
3133	cIM cf. FLott END
3134	cIM
3135	IF (ip_ebil_phy.ge.2) THEN
3136	ztit='after orography'
3137	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3138	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3139	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3140	END IF
3141	c
3142	c
3143	cAA
3144	cAA Installation de l'interface online-offline pour traceurs
3145	cAA
3146	c====================================================================
3147	c Calcul des tendances traceurs
3148	c====================================================================
3149	C
3150	IF (config_inca /= 'none') rnpb=.FALSE.
3151
3152	call phytrac ( rnpb,
3153	I itap,
3154	I julien,
3155	I gmtime,
3156	I debut,
3157	I lafin,
3158	I nqmax-2,
3159	I nlon,
3160	I nlev,
3161	I dtime,
3162	I u,
3163	I v,
3164	I t,
3165	I paprs,
3166	I pplay,
3167	I pmfu,
3168	I pmfd,
3169	I pen_u,
3170	I pde_u,
3171	I pen_d,
3172	I pde_d,
3173	I ycoefh,
3174	I fm_therm,
3175	I entr_therm,
3176	I yu1,
3177	I yv1,
3178	I ftsol,
3179	I pctsrf,
3180	I rlat,
3181	I frac_impa,
3182	I frac_nucl,
3183	I rlon,
3184	I presnivs,
3185	I pphis,
3186	I pphi,
3187	I albsol1,
3188	I qx(1,1,1),
3189	I rhcl,
3190	I cldfra,
3191	I rneb,
3192	I diafra,
3193	I cldliq,
3194	I itop_con,
3195	I ibas_con,
3196	I pmflxr,
3197	I pmflxs,
3198	I prfl,
3199	I psfl,
3200	I da,
3201	I phi,
3202	I mp,
3203	I upwd,
3204	I dnwd,
3205	I aerosol_couple,
3206	I flxmass_w,
3207	I tau_inca,
3208	I piz_inca,
3209	I cg_inca,
3210	I ccm,
3211	I rfname,
3212	O tr_seri)
3213
3214	IF (offline) THEN
3215
3216	print*,'Attention on met a 0 les thermiques pour phystoke'
3217	call phystokenc (
3218	I nlon,nlev,pdtphys,rlon,rlat,
3219	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3220	I fm_therm,entr_therm,
3221	I ycoefh,yu1,yv1,ftsol,pctsrf,
3222	I frac_impa, frac_nucl,
3223	I pphis,airephy,dtime,itap)
3224
3225
3226	ENDIF
3227
3228	c
3229	c Calculer le transport de l'eau et de l'energie (diagnostique)
3230	c
3231	CALL transp (paprs,zxtsol,
3232	e t_seri, q_seri, u_seri, v_seri, zphi,
3233	s ve, vq, ue, uq)
3234	c
3235	cIM global posePB BEG
3236	IF(1.EQ.0) THEN
3237	c
3238	CALL transp_lay (paprs,zxtsol,
3239	e t_seri, q_seri, u_seri, v_seri, zphi,
3240	s ve_lay, vq_lay, ue_lay, uq_lay)
3241	c
3242	ENDIF !(1.EQ.0) THEN
3243	cIM global posePB END
3244	c Accumuler les variables a stocker dans les fichiers histoire:
3245	c
3246	c+jld ec_conser
3247	DO k = 1, klev
3248	DO i = 1, klon
3249	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3250	d_t_ec(i,k)=0.5/ZRCPD
3251	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3252	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3253	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3254	END DO
3255	END DO
3256	c-jld ec_conser
3257	cIM
3258	IF (ip_ebil_phy.ge.1) THEN
3259	ztit='after physic'
3260	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3261	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3262	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3263	C Comme les tendances de la physique sont ajoute dans la dynamique,
3264	C on devrait avoir que la variation d'entalpie par la dynamique
3265	C est egale a la variation de la physique au pas de temps precedent.
3266	C Donc la somme de ces 2 variations devrait etre nulle.
3267	call diagphy(airephy,ztit,ip_ebil_phy
3268	e , topsw, toplw, solsw, sollw, sens
3269	e , evap, rain_fall, snow_fall, ztsol
3270	e , d_h_vcol, d_qt, d_ec
3271	s , fs_bound, fq_bound )
3272	C
3273	d_h_vcol_phy=d_h_vcol
3274	C
3275	END IF
3276	C
3277	c=======================================================================
3278	c SORTIES
3279	c=======================================================================
3280
3281	cIM Interpolation sur les niveaux de pression du NMC
3282	c -------------------------------------------------
3283	c
3284	#include "calcul_STDlev.h"
3285	c
3286	c slp sea level pressure
3287	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3288	c
3289	ccc prw = eau precipitable
3290	DO i = 1, klon
3291	prw(i) = 0.
3292	DO k = 1, klev
3293	prw(i) = prw(i) +
3294	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3295	ENDDO
3296	ENDDO
3297	c
3298	cIM initialisation + calculs divers diag AMIP2
3299	c
3300	#include "calcul_divers.h"
3301	c
3302	IF (config_inca /= 'none') THEN
3303	#ifdef INCA
3304	CALL VTe(VTphysiq)
3305	CALL VTb(VTinca)
3306
3307	CALL chemhook_end (calday,
3308	$ dtime,
3309	$ pplay,
3310	$ t_seri,
3311	$ tr_seri,
3312	$ nbtr,
3313	$ paprs,
3314	$ q_seri,
3315	$ annee_ref,
3316	$ day_ini,
3317	$ airephy,
3318	$ xjour,
3319	$ pphi,
3320	$ pphis,
3321	$ zx_rh)
3322
3323	CALL VTe(VTinca)
3324	CALL VTb(VTphysiq)
3325	#endif
3326	END IF
3327
3328	c=============================================================
3329	c
3330	c Convertir les incrementations en tendances
3331	c
3332	if (mydebug) then
3333	call writefield_phy('u_seri',u_seri,llm)
3334	call writefield_phy('v_seri',v_seri,llm)
3335	call writefield_phy('t_seri',t_seri,llm)
3336	call writefield_phy('q_seri',q_seri,llm)
3337	endif
3338
3339	DO k = 1, klev
3340	DO i = 1, klon
3341	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3342	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3343	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3344	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3345	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3346	ENDDO
3347	ENDDO
3348	c
3349	IF (nqmax.GE.3) THEN
3350	DO iq = 3, nqmax
3351	DO k = 1, klev
3352	DO i = 1, klon
3353	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3354	ENDDO
3355	ENDDO
3356	ENDDO
3357	ENDIF
3358	c
3359	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3360	cIM global posePB#include "write_bilKP_ins.h"
3361	cIM global posePB#include "write_bilKP_ave.h"
3362	c
3363	c Sauvegarder les valeurs de t et q a la fin de la physique:
3364	c
3365	DO k = 1, klev
3366	DO i = 1, klon
3367	t_ancien(i,k) = t_seri(i,k)
3368	q_ancien(i,k) = q_seri(i,k)
3369	ENDDO
3370	ENDDO
3371	c
3372	!==========================================================================
3373	! Sorties des tendances pour un point particulier
3374	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3375	! pour le debug
3376	! La valeur de igout est attribuee plus haut dans le programme
3377	!==========================================================================
3378
3379	if (prt_level.ge.1) then
3380	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3381	write(lunout,*)
3382	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
3383	write(lunout,*)
3384	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys,
3385	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3386	s pctsrf(igout,is_sic)
3387	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3388	do k=1,nlev
3389	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3390	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3391	s d_t_eva(igout,k)
3392	enddo
3393	write(lunout,*) 'cool,heat'
3394	do k=1,nlev
3395	write(lunout,*) cool(igout,k),heat(igout,k)
3396	enddo
3397
3398	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3399	do k=1,nlev
3400	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3401	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3402	enddo
3403
3404	write(lunout,*) 'd_ps ',d_ps(igout)
3405	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3406	do k=1,nlev
3407	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3408	s d_qx(igout,k,1),d_qx(igout,k,2)
3409	enddo
3410	endif
3411
3412	!==========================================================================
3413
3414	c============================================================
3415	c Calcul de la temperature potentielle
3416	c============================================================
3417	DO k = 1, klev
3418	DO i = 1, klon
3419	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3420	ENDDO
3421	ENDDO
3422	c
3423
3424	c 22.03.04 BEG
3425	c=============================================================
3426	c Ecriture des sorties
3427	c=============================================================
3428	#ifdef CPP_IOIPSL
3429
3430	c Recupere des varibles calcule dans differents modules
3431	c pour ecriture dans histxxx.nc
3432
3433	! Get some variables from module fonte_neige_mod
3434	CALL fonte_neige_get_vars(pctsrf,
3435	. zxfqcalving, zxfqfonte, zxffonte)
3436
3437
3438	c Commente par abderrahmane le 11 2 08
3439	c#ifdef histhf
3440	c#include "write_histhf.h"
3441	c#endif
3442
3443	c#ifdef histday
3444	c#include "write_histday.h"
3445	c#endif
3446
3447	c#ifdef histmth
3448	c#include "write_histmth.h"
3449	c#endif
3450
3451	c#ifdef histins
3452	c#include "write_histins.h"
3453	c#endif
3454
3455	#include "phys_output_write.h"
3456
3457	#ifdef histISCCP
3458	#include "write_histISCCP.h"
3459	#endif
3460
3461	#ifdef histmthNMC
3462	#include "write_histmthNMC.h"
3463	#endif
3464
3465	#include "write_histday_seri.h"
3466
3467	#include "write_paramLMDZ_phy.h"
3468
3469	#endif
3470
3471	c 22.03.04 END
3472	c
3473	c====================================================================
3474	c Si c'est la fin, il faut conserver l'etat de redemarrage
3475	c====================================================================
3476	c
3477
3478
3479	IF (lafin) THEN
3480	itau_phy = itau_phy + itap
3481	CALL phyredem ("restartphy.nc")
3482	! open(97,form="unformatted",file="finbin")
3483	! write(97) u_seri,v_seri,t_seri,q_seri
3484	! close(97)
3485	ENDIF
3486
3487
3488	RETURN
3489	END
3490	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3491	IMPLICIT none
3492	c
3493	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3494	c la conservation de l'eau
3495	c
3496	#include "YOMCST.h"
3497	INTEGER klon,klev
3498	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3499	REAL aire(klon)
3500	REAL qtotal, zx, qcheck
3501	INTEGER i, k
3502	c
3503	zx = 0.0
3504	DO i = 1, klon
3505	zx = zx + aire(i)
3506	ENDDO
3507	qtotal = 0.0
3508	DO k = 1, klev
3509	DO i = 1, klon
3510	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3511	. *(paprs(i,k)-paprs(i,k+1))/RG
3512	ENDDO
3513	ENDDO
3514	c
3515	qcheck = qtotal/zx
3516	c
3517	RETURN
3518	END
3519	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3520	IMPLICIT none
3521	c
3522	c Tranformer une variable de la grille physique a
3523	c la grille d'ecriture
3524	c
3525	INTEGER nfield,nlon,iim,jjmp1, jjm
3526	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3527	c
3528	INTEGER i, n, ig
3529	c
3530	jjm = jjmp1 - 1
3531	DO n = 1, nfield
3532	DO i=1,iim
3533	ecrit(i,n) = fi(1,n)
3534	ecrit(i+jjm*iim,n) = fi(nlon,n)
3535	ENDDO
3536	DO ig = 1, nlon - 2
3537	ecrit(iim+ig,n) = fi(1+ig,n)
3538	ENDDO
3539	ENDDO
3540	RETURN
3541	END

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

Download in other formats:

Original Format