Context Navigation

physiq.F @ 1422

Last change on this file since 1422 was 1422, checked in by jghattas, 14 years ago
Corrected error in use of flag prt_level. MPL
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 122.7 KB

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

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

Download in other formats:

Original Format