Context Navigation

physiq.F @ 1604

Last change on this file since 1604 was 1604, checked in by lguez, 12 years ago

Removed two "fi" with no corresponding "if" in "makdim".

In procedure "coefkz", when calculating cloud fraction "zfr", "zq" can
be zero at high altitude. Avoid division by zero, set "zfr" to 0 when
"zq" is 0.

In "physiq", allocatable arrays "tabijgcm", "longcm", "latgcm",
"igcm", "jgcm" cannot be arguments of "phys_output_open" if they have
not been allocated. Allocate them with zero size when
'npCFMIP_param.data' cannot be opened.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 137.0 KB

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

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

Download in other formats:

Original Format