Context Navigation

physiq.F @ 1663

Last change on this file since 1663 was 1641, checked in by jghattas, 12 years ago

Corrected bug while using parameter phys_out_filetimestep instead of ecrit_mth : cosp did not write at correct timestep. This correction also corrects timestep in histday_seri.nc file.

Replaced use of lev_histmth by logical swaero_diag in sw_aeroAR4.f90. swaero_diag is true if at least one of following variables are defined in output file : topswad, topswai, solswad, solswai

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

File size: 141.7 KB

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

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

Download in other formats:

Original Format