Context Navigation

physiq.F @ 1613

Last change on this file since 1613 was 1613, checked in by lguez, 13 years ago
Some input arguments of aaam_bud are defined only if ok_orodr.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 136.8 KB

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

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

Download in other formats:

Original Format