Context Navigation

physiq.F @ 1527

Last change on this file since 1527 was 1527, checked in by idelkadi, 14 years ago
Correction d'erreur introduite dans la version 1526
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 129.2 KB

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

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

Download in other formats:

Original Format