Context Navigation

physiq.F @ 1411

Last change on this file since 1411 was 1411, checked in by idelkadi, 14 years ago
Nettoyage dans physiq.F des anciennes versions versions des options iflag_cldcon=5 ou 6. Dans isrtilp.F, iflag_cldcon=5 : la version bi-gaussienne partout et iflag_cldcon=6 : bi-gaussiennes pour les thermiques et la lognormale ailleurs
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 122.5 KB

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

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

Download in other formats:

Original Format