Context Navigation

physiq.F @ 1538

Last change on this file since 1538 was 1538, checked in by musat, 13 years ago

Ajout fichier 1D mensuel paramLMDZ_phy.nc contenant

les parametres orbitaux
les taux des GES
les moyennes globales (bils, evap, evap_land, flat, nettop0, nettop, precip, tsol, t2m, prw)

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

File size: 130.0 KB

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

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

Download in other formats:

Original Format