Context Navigation

physiq.F @ 1563

Last change on this file since 1563 was 1563, checked in by jghattas, 13 years ago

Added reading of paramter type_trac from *.def file : type_trac=lmdz(default), inca or repr(soon). While running with INCA, 2 parameters are now necessare in .def : type_trac=inca and config_inca=aero/chem. If type_trac=lmdz or repr, config_inca will not be used.

Removed print of ecrit_mth, ecrit_day etc in physiq.F. Removed the variable ecrit_hf2mth which is no longer used.

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

File size: 136.0 KB

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

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

Download in other formats:

Original Format