Context Navigation

physiq.F @ 1516

Last change on this file since 1516 was 1516, checked in by idelkadi, 13 years ago

Modification dans physiq.F pour que les pluies stratiformes n'alimentent les pochers que si le spoches occupent plus de 10% de la surface.
Dans le cas iflag_wake=2, on ne prend en compte les tendances de fisrtilp pour forcer les poches que si les poches existent deja (en controlant si elles couvrent deja 10% de la maille)
On autorise differentes formulations pour la vitesse a la base de la convection WB (plus precisement au niveau de convection libre LFC, en haut de la zone d'inhibition) en fonction de flag_wb, lu dans conv_param.data

flag_wb=0 : wb=wbmax, wbmax etant lu dans conv_param.data
flag_wb=1 : wb=wbmax/[1+ 500hPa / (ps-p_LFC) ]
Concerne cv3_routines.F cv3param.h cv3p1_closure.F

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

File size: 128.8 KB

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

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

Download in other formats:

Original Format