Context Navigation

physiq.F @ 1639

Last change on this file since 1639 was 1639, checked in by jghattas, 12 years ago
Corrected last commit : lnblnk1 do not existe anymore. The model did not compile at mercure.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 141.9 KB

Line
1	! $Id: physiq.F 1639 2012-07-24 05:17:23Z jghattas $
2	c#define IO_DEBUG
3
4	SUBROUTINE physiq (nlon,nlev,
5	. debut,lafin,jD_cur, jH_cur,pdtphys,
6	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
7	. u,v,t,qx,
8	. flxmass_w,
9	. d_u, d_v, d_t, d_qx, d_ps
10	. , dudyn
11	. , PVteta)
12
13	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
14	$ histwrite, ju2ymds, ymds2ju, ioget_year_len
15	USE comgeomphy
16	USE phys_cal_mod
17	USE write_field_phy
18	USE dimphy
19	USE infotrac
20	USE mod_grid_phy_lmdz
21	USE mod_phys_lmdz_para
22	USE iophy
23	USE misc_mod, mydebug=>debug
24	USE vampir
25	USE pbl_surface_mod, ONLY : pbl_surface
26	USE change_srf_frac_mod
27	USE surface_data, ONLY : type_ocean, ok_veget
28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
29	USE phys_state_var_mod ! Variables sauvegardees de la physique
30	USE phys_output_var_mod ! Variables pour les ecritures des sorties
31	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
32	USE phys_output_mod
33	use open_climoz_m, only: open_climoz ! ozone climatology from a file
34	use regr_pr_av_m, only: regr_pr_av
35	use netcdf95, only: nf95_close
36	cIM for NMC files
37	use netcdf, only: nf90_fill_real
38	use mod_phys_lmdz_mpi_data, only: is_mpi_root
39	USE aero_mod
40	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
41	use conf_phys_m, only: conf_phys
42	use radlwsw_m, only: radlwsw
43	USE control_mod
44	#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
1664
1665	ecrit_hf = ecrit_hf * un_jour
1666	cIM
1667	IF(ecrit_day.LE.1.) THEN
1668	ecrit_day = ecrit_day * un_jour !en secondes
1669	ENDIF
1670	cIM
1671	ecrit_mth = ecrit_mth * un_jour
1672	ecrit_ins = ecrit_ins * un_jour
1673	ecrit_reg = ecrit_reg * un_jour
1674	ecrit_tra = ecrit_tra * un_jour
1675	ecrit_LES = ecrit_LES * un_jour
1676	c
1677
1678	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1679	date0 = jD_ref
1680	WRITE(,) 'physiq date0 : ',date0
1681	c
1682	c
1683	c
1684	c Prescrire l'ozone dans l'atmosphere
1685	c
1686	c
1687	cc DO i = 1, klon
1688	cc DO k = 1, klev
1689	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1690	cc ENDDO
1691	cc ENDDO
1692	c
1693	IF (type_trac == 'inca') THEN
1694	#ifdef INCA
1695	CALL VTe(VTphysiq)
1696	CALL VTb(VTinca)
1697	! iii = MOD(NINT(xjour),360)
1698	! calday = REAL(iii) + jH_cur
1699	calday = REAL(days_elapsed) + jH_cur
1700	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1701
1702	CALL chemini(
1703	$ rg,
1704	$ ra,
1705	$ airephy,
1706	$ rlat,
1707	$ rlon,
1708	$ presnivs,
1709	$ calday,
1710	$ klon,
1711	$ nqtot,
1712	$ pdtphys,
1713	$ annee_ref,
1714	$ day_ref,
1715	$ itau_phy)
1716
1717	CALL VTe(VTinca)
1718	CALL VTb(VTphysiq)
1719	#endif
1720	END IF
1721	c
1722	c
1723	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1724	! Nouvelle initialisation pour le rayonnement RRTM
1725	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1726
1727	call iniradia(klon,klev,paprs(1,1:klev+1))
1728
1729	C$omp single
1730	if (read_climoz >= 1) then
1731	call open_climoz(ncid_climoz, press_climoz)
1732	END IF
1733	C$omp end single
1734	c
1735	cIM betaCRF
1736	pfree=70000. !Pa
1737	beta_pbl=1.
1738	beta_free=1.
1739	lon1_beta=-180.
1740	lon2_beta=+180.
1741	lat1_beta=90.
1742	lat2_beta=-90.
1743	mskocean_beta=.FALSE.
1744
1745	OPEN(99,file='beta_crf.data',status='old',
1746	$ form='formatted',err=9999)
1747	READ(99,*,end=9998) pfree
1748	READ(99,*,end=9998) beta_pbl
1749	READ(99,*,end=9998) beta_free
1750	READ(99,*,end=9998) lon1_beta
1751	READ(99,*,end=9998) lon2_beta
1752	READ(99,*,end=9998) lat1_beta
1753	READ(99,*,end=9998) lat2_beta
1754	READ(99,*,end=9998) mskocean_beta
1755	9998 Continue
1756	CLOSE(99)
1757	9999 Continue
1758	WRITE(,)'pfree=',pfree
1759	WRITE(,)'beta_pbl=',beta_pbl
1760	WRITE(,)'beta_free=',beta_free
1761	WRITE(,)'lon1_beta=',lon1_beta
1762	WRITE(,)'lon2_beta=',lon2_beta
1763	WRITE(,)'lat1_beta=',lat1_beta
1764	WRITE(,)'lat2_beta=',lat2_beta
1765	WRITE(,)'mskocean_beta=',mskocean_beta
1766	ENDIF
1767	!
1768	! ************** Fin de IF ( debut ) *************
1769	!
1770	!
1771	! Incrementer le compteur de la physique
1772	!
1773	itap = itap + 1
1774	c
1775	!
1776	! Update fraction of the sub-surfaces (pctsrf) and
1777	! initialize, where a new fraction has appeared, all variables depending
1778	! on the surface fraction.
1779	!
1780	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1781	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1782
1783
1784	! Update time and other variables in Reprobus
1785	IF (type_trac == 'repr') THEN
1786	#ifdef REPROBUS
1787	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1788	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1789	CALL Rtime(debut)
1790	#endif
1791	END IF
1792
1793
1794	! Tendances bidons pour les processus qui n'affectent pas certaines
1795	! variables.
1796	du0(:,:)=0.
1797	dv0(:,:)=0.
1798	dq0(:,:)=0.
1799	dql0(:,:)=0.
1800	c
1801	c Mettre a zero des variables de sortie (pour securite)
1802	c
1803	DO i = 1, klon
1804	d_ps(i) = 0.0
1805	ENDDO
1806	DO k = 1, klev
1807	DO i = 1, klon
1808	d_t(i,k) = 0.0
1809	d_u(i,k) = 0.0
1810	d_v(i,k) = 0.0
1811	ENDDO
1812	ENDDO
1813	DO iq = 1, nqtot
1814	DO k = 1, klev
1815	DO i = 1, klon
1816	d_qx(i,k,iq) = 0.0
1817	ENDDO
1818	ENDDO
1819	ENDDO
1820	da(:,:)=0.
1821	mp(:,:)=0.
1822	phi(:,:,:)=0.
1823	c
1824	c Ne pas affecter les valeurs entrees de u, v, h, et q
1825	c
1826	DO k = 1, klev
1827	DO i = 1, klon
1828	t_seri(i,k) = t(i,k)
1829	u_seri(i,k) = u(i,k)
1830	v_seri(i,k) = v(i,k)
1831	q_seri(i,k) = qx(i,k,ivap)
1832	ql_seri(i,k) = qx(i,k,iliq)
1833	qs_seri(i,k) = 0.
1834	ENDDO
1835	ENDDO
1836	IF (nqtot.GE.3) THEN
1837	DO iq = 3, nqtot
1838	DO k = 1, klev
1839	DO i = 1, klon
1840	tr_seri(i,k,iq-2) = qx(i,k,iq)
1841	ENDDO
1842	ENDDO
1843	ENDDO
1844	ELSE
1845	DO k = 1, klev
1846	DO i = 1, klon
1847	tr_seri(i,k,1) = 0.0
1848	ENDDO
1849	ENDDO
1850	ENDIF
1851	C
1852	DO i = 1, klon
1853	ztsol(i) = 0.
1854	ENDDO
1855	DO nsrf = 1, nbsrf
1856	DO i = 1, klon
1857	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1858	ENDDO
1859	ENDDO
1860	cIM
1861	IF (ip_ebil_phy.ge.1) THEN
1862	ztit='after dynamic'
1863	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1864	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1865	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1866	C Comme les tendances de la physique sont ajoute dans la dynamique,
1867	C on devrait avoir que la variation d'entalpie par la dynamique
1868	C est egale a la variation de la physique au pas de temps precedent.
1869	C Donc la somme de ces 2 variations devrait etre nulle.
1870	call diagphy(airephy,ztit,ip_ebil_phy
1871	e , zero_v, zero_v, zero_v, zero_v, zero_v
1872	e , zero_v, zero_v, zero_v, ztsol
1873	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1874	s , fs_bound, fq_bound )
1875	END IF
1876
1877	c Diagnostiquer la tendance dynamique
1878	c
1879	IF (ancien_ok) THEN
1880	DO k = 1, klev
1881	DO i = 1, klon
1882	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1883	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1884	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1885	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1886	ENDDO
1887	ENDDO
1888	ELSE
1889	DO k = 1, klev
1890	DO i = 1, klon
1891	d_u_dyn(i,k) = 0.0
1892	d_v_dyn(i,k) = 0.0
1893	d_t_dyn(i,k) = 0.0
1894	d_q_dyn(i,k) = 0.0
1895	ENDDO
1896	ENDDO
1897	ancien_ok = .TRUE.
1898	ENDIF
1899	c
1900	c Ajouter le geopotentiel du sol:
1901	c
1902	DO k = 1, klev
1903	DO i = 1, klon
1904	zphi(i,k) = pphi(i,k) + pphis(i)
1905	ENDDO
1906	ENDDO
1907	c
1908	c Verifier les temperatures
1909	c
1910	cIM BEG
1911	IF (check) THEN
1912	amn=MIN(ftsol(1,is_ter),1000.)
1913	amx=MAX(ftsol(1,is_ter),-1000.)
1914	DO i=2, klon
1915	amn=MIN(ftsol(i,is_ter),amn)
1916	amx=MAX(ftsol(i,is_ter),amx)
1917	ENDDO
1918	c
1919	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1920	ENDIF !(check) THEN
1921	cIM END
1922	c
1923	CALL hgardfou(t_seri,ftsol,'debutphy')
1924	c
1925	cIM BEG
1926	IF (check) THEN
1927	amn=MIN(ftsol(1,is_ter),1000.)
1928	amx=MAX(ftsol(1,is_ter),-1000.)
1929	DO i=2, klon
1930	amn=MIN(ftsol(i,is_ter),amn)
1931	amx=MAX(ftsol(i,is_ter),amx)
1932	ENDDO
1933	c
1934	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1935	ENDIF !(check) THEN
1936	cIM END
1937	c
1938	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1939	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1940	c
1941	if (read_climoz >= 1) then
1942	C Ozone from a file
1943	! Update required ozone index:
1944	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1945	$ / ioget_year_len(year_cur) * 360.) + 1
1946	if (ro3i == 361) ro3i = 360
1947	C (This should never occur, except perhaps because of roundup
1948	C error. See documentation.)
1949	if (ro3i /= co3i) then
1950	C Update ozone field:
1951	if (read_climoz == 1) then
1952	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1953	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1954	else
1955	C read_climoz == 2
1956	call regr_pr_av(ncid_climoz,
1957	$ (/"tro3 ", "tro3_daylight"/),
1958	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1959	$ v3=wo)
1960	end if
1961	! Convert from mole fraction of ozone to column density of ozone in a
1962	! cell, in kDU:
1963	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1964	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1965	C (By regridding ozone values for LMDZ only once every 360th of
1966	C year, we have already neglected the variation of pressure in one
1967	C 360th of year. So do not recompute "wo" at each time step even if
1968	C "zmasse" changes a little.)
1969	co3i = ro3i
1970	end if
1971	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1972	C Once per day, update ozone from Royer:
1973	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1974	ENDIF
1975	c
1976	c Re-evaporer l'eau liquide nuageuse
1977	c
1978	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1979	DO i = 1, klon
1980	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1981	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1982	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1983	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1984	zb = MAX(0.0,ql_seri(i,k))
1985	za = - MAX(0.0,ql_seri(i,k))
1986	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1987	t_seri(i,k) = t_seri(i,k) + za
1988	q_seri(i,k) = q_seri(i,k) + zb
1989	ql_seri(i,k) = 0.0
1990	d_t_eva(i,k) = za
1991	d_q_eva(i,k) = zb
1992	ENDDO
1993	ENDDO
1994	cIM
1995	IF (ip_ebil_phy.ge.2) THEN
1996	ztit='after reevap'
1997	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1998	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1999	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2000	call diagphy(airephy,ztit,ip_ebil_phy
2001	e , zero_v, zero_v, zero_v, zero_v, zero_v
2002	e , zero_v, zero_v, zero_v, ztsol
2003	e , d_h_vcol, d_qt, d_ec
2004	s , fs_bound, fq_bound )
2005	C
2006	END IF
2007
2008	c
2009	c=========================================================================
2010	! Calculs de l'orbite.
2011	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2012	! doit donc etre placé avant radlwsw et pbl_surface
2013
2014	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2015	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2016	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2017	!
2018	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2019	! solarlong0
2020	if (solarlong0<-999.) then
2021	if (new_orbit) then
2022	! calcul selon la routine utilisee pour les planetes
2023	call solarlong(day_since_equinox, zlongi, dist)
2024	else
2025	! calcul selon la routine utilisee pour l'AR4
2026	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2027	endif
2028	else
2029	zlongi=solarlong0 ! longitude solaire vraie
2030	dist=1. ! distance au soleil / moyenne
2031	endif
2032	if(prt_level.ge.1) &
2033	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2034
2035
2036	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2037	! Calcul de l'ensoleillement :
2038	! ============================
2039	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2040	! l'annee a partir d'une formule analytique.
2041	! Cet ensoleillement est symmétrique autour de l'équateur et
2042	! non nul aux poles.
2043	IF (abs(solarlong0-1000.)<1.e-4) then
2044	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
2045	ELSE
2046	! Avec ou sans cycle diurne
2047	IF (cycle_diurne) THEN
2048	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
2049	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
2050	ELSE
2051	CALL angle(zlongi, rlat, fract, rmu0)
2052	ENDIF
2053	ENDIF
2054
2055	if (mydebug) then
2056	call writefield_phy('u_seri',u_seri,llm)
2057	call writefield_phy('v_seri',v_seri,llm)
2058	call writefield_phy('t_seri',t_seri,llm)
2059	call writefield_phy('q_seri',q_seri,llm)
2060	endif
2061
2062	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2063	c Appel au pbl_surface : Planetary Boudary Layer et Surface
2064	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2065	c turbulent du couche limit.
2066	c
2067	c Certains varibales de sorties de pbl_surface sont utiliser que pour
2068	c ecriture des fihiers hist_XXXX.nc, ces sont :
2069	c qsol, zq2m, s_pblh, s_lcl,
2070	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2071	c s_therm, s_trmb1, s_trmb2, s_trmb3,
2072	c zxrugs, zu10m, zv10m, fder,
2073	c zxqsurf, rh2m, zxfluxu, zxfluxv,
2074	c frugs, agesno, fsollw, fsolsw,
2075	c d_ts, fevap, fluxlat, t2m,
2076	c wfbils, wfbilo, fluxt, fluxu, fluxv,
2077	c
2078	c Certains ne sont pas utiliser du tout :
2079	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2080	c
2081
2082	if (iflag_pbl/=0) then
2083
2084	CALL pbl_surface(
2085	e dtime, date0, itap, days_elapsed+1,
2086	e debut, lafin,
2087	e rlon, rlat, rugoro, rmu0,
2088	e rain_fall, snow_fall, solsw, sollw,
2089	e t_seri, q_seri, u_seri, v_seri,
2090	e pplay, paprs, pctsrf,
2091	+ ftsol, falb1, falb2, u10m, v10m,
2092	s sollwdown, cdragh, cdragm, u1, v1,
2093	s albsol1, albsol2, sens, evap,
2094	s zxtsol, zxfluxlat, zt2m, qsat2m,
2095	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
2096	s coefh, coefm, slab_wfbils,
2097	d qsol, zq2m, s_pblh, s_lcl,
2098	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2099	d s_therm, s_trmb1, s_trmb2, s_trmb3,
2100	d zxrugs, zu10m, zv10m, fder,
2101	d zxqsurf, rh2m, zxfluxu, zxfluxv,
2102	d frugs, agesno, fsollw, fsolsw,
2103	d d_ts, fevap, fluxlat, t2m,
2104	d wfbils, wfbilo, fluxt, fluxu, fluxv,
2105	- dsens, devap, zxsnow,
2106	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
2107
2108
2109	!-----------------------------------------------------------------------------------------
2110	! ajout des tendances de la diffusion turbulente
2111	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
2112	!-----------------------------------------------------------------------------------------
2113
2114	if (mydebug) then
2115	call writefield_phy('u_seri',u_seri,llm)
2116	call writefield_phy('v_seri',v_seri,llm)
2117	call writefield_phy('t_seri',t_seri,llm)
2118	call writefield_phy('q_seri',q_seri,llm)
2119	endif
2120
2121
2122	IF (ip_ebil_phy.ge.2) THEN
2123	ztit='after surface_main'
2124	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2125	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2126	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2127	call diagphy(airephy,ztit,ip_ebil_phy
2128	e , zero_v, zero_v, zero_v, zero_v, sens
2129	e , evap , zero_v, zero_v, ztsol
2130	e , d_h_vcol, d_qt, d_ec
2131	s , fs_bound, fq_bound )
2132	END IF
2133
2134	ENDIF
2135	c =================================================================== c
2136	c Calcul de Qsat
2137
2138	DO k = 1, klev
2139	DO i = 1, klon
2140	zx_t = t_seri(i,k)
2141	IF (thermcep) THEN
2142	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2143	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2144	zx_qs = MIN(0.5,zx_qs)
2145	zcor = 1./(1.-retv*zx_qs)
2146	zx_qs = zx_qs*zcor
2147	ELSE
2148	IF (zx_t.LT.t_coup) THEN
2149	zx_qs = qsats(zx_t)/pplay(i,k)
2150	ELSE
2151	zx_qs = qsatl(zx_t)/pplay(i,k)
2152	ENDIF
2153	ENDIF
2154	zqsat(i,k)=zx_qs
2155	ENDDO
2156	ENDDO
2157
2158	if (prt_level.ge.1) then
2159	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2160	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2161	endif
2162	c
2163	c Appeler la convection (au choix)
2164	c
2165	DO k = 1, klev
2166	DO i = 1, klon
2167	conv_q(i,k) = d_q_dyn(i,k)
2168	. + d_q_vdf(i,k)/dtime
2169	conv_t(i,k) = d_t_dyn(i,k)
2170	. + d_t_vdf(i,k)/dtime
2171	ENDDO
2172	ENDDO
2173	IF (check) THEN
2174	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2175	WRITE(lunout,*) "avantcon=", za
2176	ENDIF
2177	zx_ajustq = .FALSE.
2178	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2179	IF (zx_ajustq) THEN
2180	DO i = 1, klon
2181	z_avant(i) = 0.0
2182	ENDDO
2183	DO k = 1, klev
2184	DO i = 1, klon
2185	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2186	. *(paprs(i,k)-paprs(i,k+1))/RG
2187	ENDDO
2188	ENDDO
2189	ENDIF
2190
2191	c Calcule de vitesse verticale a partir de flux de masse verticale
2192	DO k = 1, klev
2193	DO i = 1, klon
2194	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2195	END DO
2196	END DO
2197	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
2198	$ omega(igout, :)
2199
2200	IF (iflag_con.EQ.1) THEN
2201	abort_message ='reactiver le call conlmd dans physiq.F'
2202	CALL abort_gcm (modname,abort_message,1)
2203	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2204	c . d_t_con, d_q_con,
2205	c . rain_con, snow_con, ibas_con, itop_con)
2206	ELSE IF (iflag_con.EQ.2) THEN
2207	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2208	e conv_t, conv_q, -evap, omega,
2209	s d_t_con, d_q_con, rain_con, snow_con,
2210	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2211	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2212	d_u_con = 0.
2213	d_v_con = 0.
2214
2215	WHERE (rain_con < 0.) rain_con = 0.
2216	WHERE (snow_con < 0.) snow_con = 0.
2217	DO i = 1, klon
2218	ibas_con(i) = klev+1 - kcbot(i)
2219	itop_con(i) = klev+1 - kctop(i)
2220	ENDDO
2221	ELSE IF (iflag_con.GE.3) THEN
2222	c nb of tracers for the KE convection:
2223	c MAF la partie traceurs est faite dans phytrac
2224	c on met ntra=1 pour limiter les appels mais on peut
2225	c supprimer les calculs / ftra.
2226	ntra = 1
2227
2228	c=====================================================================================
2229	cajout pour la parametrisation des poches froides:
2230	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2231	do k=1,klev
2232	do i=1,klon
2233	if (iflag_wake>=1) then
2234	t_wake(i,k) = t_seri(i,k)
2235	. +(1-wake_s(i))*wake_deltat(i,k)
2236	q_wake(i,k) = q_seri(i,k)
2237	. +(1-wake_s(i))*wake_deltaq(i,k)
2238	t_undi(i,k) = t_seri(i,k)
2239	. -wake_s(i)*wake_deltat(i,k)
2240	q_undi(i,k) = q_seri(i,k)
2241	. -wake_s(i)*wake_deltaq(i,k)
2242	else
2243	t_wake(i,k) = t_seri(i,k)
2244	q_wake(i,k) = q_seri(i,k)
2245	t_undi(i,k) = t_seri(i,k)
2246	q_undi(i,k) = q_seri(i,k)
2247	endif
2248	enddo
2249	enddo
2250
2251	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2252	cc-- pour le soulevement des particules dans le modele convectif
2253	c
2254	do i = 1,klon
2255	ALE(i) = 0.
2256	ALP(i) = 0.
2257	enddo
2258	c
2259	ccalcul de ale_wake et alp_wake
2260	if (iflag_wake>=1) then
2261	if (itap .le. it_wape_prescr) then
2262	do i = 1,klon
2263	ale_wake(i) = wape_prescr
2264	alp_wake(i) = fip_prescr
2265	enddo
2266	else
2267	do i = 1,klon
2268	cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2269	ccc ale_wake(i) = 0.5wake_cstar(i)*2
2270	ale_wake(i) = wake_pe(i)
2271	alp_wake(i) = wake_fip(i)
2272	enddo
2273	endif
2274	else
2275	do i = 1,klon
2276	ale_wake(i) = 0.
2277	alp_wake(i) = 0.
2278	enddo
2279	endif
2280	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2281	cdans le thermique sinon
2282	if (iflag_coupl.eq.0) then
2283	if (debut.and.prt_level.gt.9)
2284	$ WRITE(lunout,*)'ALE et ALP imposes'
2285	do i = 1,klon
2286	con ne couple que ale
2287	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2288	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2289	con ne couple que alp
2290	c ALP(i) = alp_wake(i) + Alp_bl(i)
2291	ALP(i) = alp_wake(i) + alp_bl_prescr
2292	enddo
2293	else
2294	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2295	! do i = 1,klon
2296	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2297	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2298	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2299	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2300	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2301	! enddo
2302
2303	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2304	! Modif FH 2010/04/27. Sans doute temporaire.
2305	! Deux options pour le alp_offset : constant si >Ã 0 ou proportionnel Ãa
2306	! w si <0
2307	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2308	do i = 1,klon
2309	ALE(i) = max(ale_wake(i),Ale_bl(i))
2310	ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
2311	if (iflag_trig_bl.ge.1) then
2312	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2313	endif
2314	ccc fin nrlmd le 10/04/2012
2315	if (alp_offset>=0.) then
2316	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2317	else
2318	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2319	if (alp(i)<0.) then
2320	print*,'ALP ',alp(i),alp_wake(i)
2321	s ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2322	endif
2323	endif
2324	enddo
2325	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2326
2327	endif
2328	do i=1,klon
2329	if (alp(i)>alp_max) then
2330	IF(prt_level>9)WRITE(lunout,*) &
2331	& 'WARNING SUPER ALP (seuil=',alp_max,
2332	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2333	alp(i)=alp_max
2334	endif
2335	if (ale(i)>ale_max) then
2336	IF(prt_level>9)WRITE(lunout,*) &
2337	& 'WARNING SUPER ALE (seuil=',ale_max,
2338	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2339	ale(i)=ale_max
2340	endif
2341	enddo
2342
2343	cfin calcul ale et alp
2344	c=================================================================================================
2345
2346
2347	c sb, oct02:
2348	c Schema de convection modularise et vectorise:
2349	c (driver commun aux versions 3 et 4)
2350	c
2351	IF (ok_cvl) THEN ! new driver for convectL
2352
2353	IF (type_trac == 'repr') THEN
2354	nbtr_tmp=ntra
2355	ELSE
2356	nbtr_tmp=nbtr
2357	END IF
2358	CALL concvl (iflag_con,iflag_clos,
2359	. dtime,paprs,pplay,t_undi,q_undi,
2360	. t_wake,q_wake,wake_s,
2361	. u_seri,v_seri,tr_seri,nbtr_tmp,
2362	. ALE,ALP,
2363	. ema_work1,ema_work2,
2364	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2365	. rain_con, snow_con, ibas_con, itop_con, sigd,
2366	. ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0,
2367	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2368	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2369	. pmflxr,pmflxs,da,phi,mp,
2370	. ftd,fqd,lalim_conv,wght_th)
2371
2372	cIM begin
2373	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2374	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2375	cIM end
2376	cIM cf. FH
2377	clwcon0=qcondc
2378	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2379
2380	do i = 1, klon
2381	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2382	enddo
2383
2384	ELSE ! ok_cvl
2385
2386	c MAF conema3 ne contient pas les traceurs
2387	CALL conema3 (dtime,
2388	. paprs,pplay,t_seri,q_seri,
2389	. u_seri,v_seri,tr_seri,ntra,
2390	. ema_work1,ema_work2,
2391	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2392	. rain_con, snow_con, ibas_con, itop_con,
2393	. upwd,dnwd,dnwd0,bas,top,
2394	. Ma,cape,tvp,rflag,
2395	. pbase
2396	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2397	. ,clwcon0)
2398
2399	ENDIF ! ok_cvl
2400
2401	c
2402	c Correction precip
2403	rain_con = rain_con * cvl_corr
2404	snow_con = snow_con * cvl_corr
2405	c
2406
2407	IF (.NOT. ok_gust) THEN
2408	do i = 1, klon
2409	wd(i)=0.0
2410	enddo
2411	ENDIF
2412
2413	c =================================================================== c
2414	c Calcul des proprietes des nuages convectifs
2415	c
2416
2417	c calcul des proprietes des nuages convectifs
2418	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2419	call clouds_gno
2420	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2421
2422	c =================================================================== c
2423
2424	DO i = 1, klon
2425	itop_con(i) = min(max(itop_con(i),1),klev)
2426	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2427	ENDDO
2428
2429	DO i = 1, klon
2430	ema_pcb(i) = paprs(i,ibas_con(i))
2431	ENDDO
2432	DO i = 1, klon
2433	! L'idicage de itop_con peut cacher un pb potentiel
2434	! FH sous la dictee de JYG, CR
2435	ema_pct(i) = paprs(i,itop_con(i)+1)
2436
2437	if (itop_con(i).gt.klev-3) then
2438	if(prt_level >= 9) then
2439	write(lunout,*)'La convection monte trop haut '
2440	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2441	endif
2442	endif
2443	ENDDO
2444	ELSE IF (iflag_con.eq.0) THEN
2445	write(lunout,*) 'On n appelle pas la convection'
2446	clwcon0=0.
2447	rnebcon0=0.
2448	d_t_con=0.
2449	d_q_con=0.
2450	d_u_con=0.
2451	d_v_con=0.
2452	rain_con=0.
2453	snow_con=0.
2454	bas=1
2455	top=1
2456	ELSE
2457	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2458	CALL abort
2459	ENDIF
2460
2461	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2462	c . d_u_con, d_v_con)
2463
2464	!-----------------------------------------------------------------------------------------
2465	! ajout des tendances de la diffusion turbulente
2466	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2467	!-----------------------------------------------------------------------------------------
2468
2469	if (mydebug) then
2470	call writefield_phy('u_seri',u_seri,llm)
2471	call writefield_phy('v_seri',v_seri,llm)
2472	call writefield_phy('t_seri',t_seri,llm)
2473	call writefield_phy('q_seri',q_seri,llm)
2474	endif
2475
2476	cIM
2477	IF (ip_ebil_phy.ge.2) THEN
2478	ztit='after convect'
2479	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2480	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2481	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2482	call diagphy(airephy,ztit,ip_ebil_phy
2483	e , zero_v, zero_v, zero_v, zero_v, zero_v
2484	e , zero_v, rain_con, snow_con, ztsol
2485	e , d_h_vcol, d_qt, d_ec
2486	s , fs_bound, fq_bound )
2487	END IF
2488	C
2489	IF (check) THEN
2490	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2491	WRITE(lunout,*)"aprescon=", za
2492	zx_t = 0.0
2493	za = 0.0
2494	DO i = 1, klon
2495	za = za + airephy(i)/REAL(klon)
2496	zx_t = zx_t + (rain_con(i)+
2497	. snow_con(i))*airephy(i)/REAL(klon)
2498	ENDDO
2499	zx_t = zx_t/za*dtime
2500	WRITE(lunout,*)"Precip=", zx_t
2501	ENDIF
2502	IF (zx_ajustq) THEN
2503	DO i = 1, klon
2504	z_apres(i) = 0.0
2505	ENDDO
2506	DO k = 1, klev
2507	DO i = 1, klon
2508	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2509	. *(paprs(i,k)-paprs(i,k+1))/RG
2510	ENDDO
2511	ENDDO
2512	DO i = 1, klon
2513	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2514	. /z_apres(i)
2515	ENDDO
2516	DO k = 1, klev
2517	DO i = 1, klon
2518	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2519	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2520	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2521	ENDIF
2522	ENDDO
2523	ENDDO
2524	ENDIF
2525	zx_ajustq=.FALSE.
2526
2527	c
2528	c=============================================================================
2529	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2530	cpour la couche limite diffuse pour l instant
2531	c
2532	if (iflag_wake>=1) then
2533	DO k=1,klev
2534	DO i=1,klon
2535	dt_dwn(i,k) = ftd(i,k)
2536	wdt_PBL(i,k) = 0.
2537	dq_dwn(i,k) = fqd(i,k)
2538	wdq_PBL(i,k) = 0.
2539	M_dwn(i,k) = dnwd0(i,k)
2540	M_up(i,k) = upwd(i,k)
2541	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2542	udt_PBL(i,k) = 0.
2543	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2544	udq_PBL(i,k) = 0.
2545	ENDDO
2546	ENDDO
2547
2548	if (iflag_wake==2) then
2549	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2550	DO k = 1,klev
2551	dt_dwn(:,k)= dt_dwn(:,k)+
2552	: ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2553	dq_dwn(:,k)= dq_dwn(:,k)+
2554	: ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2555	ENDDO
2556	endif
2557	c
2558	ccalcul caracteristiques de la poche froide
2559	call calWAKE (paprs,pplay,dtime
2560	: ,t_seri,q_seri,omega
2561	: ,dt_dwn,dq_dwn,M_dwn,M_up
2562	: ,dt_a,dq_a,sigd
2563	: ,wdt_PBL,wdq_PBL
2564	: ,udt_PBL,udq_PBL
2565	o ,wake_deltat,wake_deltaq,wake_dth
2566	o ,wake_h,wake_s,wake_dens
2567	o ,wake_pe,wake_fip,wake_gfl
2568	o ,dt_wake,dq_wake
2569	o ,wake_k, t_undi,q_undi
2570	o ,wake_omgbdth,wake_dp_omgb
2571	o ,wake_dtKE,wake_dqKE
2572	o ,wake_dtPBL,wake_dqPBL
2573	o ,wake_omg,wake_dp_deltomg
2574	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2575	o ,wake_ddeltat,wake_ddeltaq)
2576	c
2577	!-----------------------------------------------------------------------------------------
2578	! ajout des tendances des poches froides
2579	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2580	! coherent avec les autres d_t_...
2581	d_t_wake(:,:)=dt_wake(:,:)*dtime
2582	d_q_wake(:,:)=dq_wake(:,:)*dtime
2583	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2584	!-----------------------------------------------------------------------------------------
2585
2586	endif
2587	c
2588	c===================================================================
2589	cJYG
2590	IF (ip_ebil_phy.ge.2) THEN
2591	ztit='after wake'
2592	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2593	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2594	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2595	call diagphy(airephy,ztit,ip_ebil_phy
2596	e , zero_v, zero_v, zero_v, zero_v, zero_v
2597	e , zero_v, zero_v, zero_v, ztsol
2598	e , d_h_vcol, d_qt, d_ec
2599	s , fs_bound, fq_bound )
2600	END IF
2601
2602	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2603	c
2604	c===================================================================
2605	c Convection seche (thermiques ou ajustement)
2606	c===================================================================
2607	c
2608	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2609	s ,seuil_inversion,weak_inversion,dthmin)
2610
2611
2612
2613	d_t_ajsb(:,:)=0.
2614	d_q_ajsb(:,:)=0.
2615	d_t_ajs(:,:)=0.
2616	d_u_ajs(:,:)=0.
2617	d_v_ajs(:,:)=0.
2618	d_q_ajs(:,:)=0.
2619	clwcon0th(:,:)=0.
2620	c
2621	c fm_therm(:,:)=0.
2622	c entr_therm(:,:)=0.
2623	c detr_therm(:,:)=0.
2624	c
2625	IF(prt_level>9)WRITE(lunout,*)
2626	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2627	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2628	if(iflag_thermals.lt.0) then
2629	c Rien
2630	c ====
2631	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2632
2633
2634	else
2635
2636	c Thermiques
2637	c ==========
2638	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2639	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2640
2641
2642	ccc nrlmd le 10/04/2012
2643	DO k=1,klev+1
2644	DO i=1,klon
2645	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2646	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2647	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2648	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2649	ENDDO
2650	ENDDO
2651	ccc fin nrlmd le 10/04/2012
2652
2653	if (iflag_thermals>=1) then
2654	call calltherm(pdtphys
2655	s ,pplay,paprs,pphi,weak_inversion
2656	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2657	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2658	s ,fm_therm,entr_therm,detr_therm
2659	s ,zqasc,clwcon0th,lmax_th,ratqscth
2660	s ,ratqsdiff,zqsatth
2661	con rajoute ale et alp, et les caracteristiques de la couche alim
2662	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
2663	s ,ztv,zpspsk,ztla,zthl
2664	ccc nrlmd le 10/04/2012
2665	e ,pbl_tke_input,pctsrf,omega,airephy
2666	s ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0
2667	s ,n2,s2,ale_bl_stat
2668	s ,therm_tke_max,env_tke_max
2669	s ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
2670	s ,alp_bl_conv,alp_bl_stat
2671	ccc fin nrlmd le 10/04/2012
2672	s )
2673
2674	ccc nrlmd le 10/04/2012
2675	c-----------Stochastic triggering-----------
2676	if (iflag_trig_bl.ge.1) then
2677	c
2678	IF (prt_level .GE. 10) THEN
2679	print *,'cin, ale_bl_stat, alp_bl_stat ',
2680	$ cin, ale_bl_stat, alp_bl_stat
2681	ENDIF
2682
2683	c----Initialisations
2684	do i=1,klon
2685	proba_notrig(i)=1.
2686	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2687	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2688	tau_trig(i)=tau_trig_shallow
2689	else
2690	tau_trig(i)=tau_trig_deep
2691	endif
2692	enddo
2693	c
2694	IF (prt_level .GE. 10) THEN
2695	print *,'random_notrig, tau_trig ',
2696	$ random_notrig, tau_trig
2697	print *,'s_trig,s2,n2 ',
2698	$ s_trig,s2,n2
2699	ENDIF
2700
2701	c----Tirage aléatoire et calcul de ale_bl_trig
2702	do i=1,klon
2703	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2704	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))**
2705	$ (n2(i)*dtime/tau_trig(i))
2706	c print *, 'proba_notrig(i) ',proba_notrig(i)
2707	if (random_notrig(i) .ge. proba_notrig(i)) then
2708	ale_bl_trig(i)=ale_bl_stat(i)
2709	else
2710	ale_bl_trig(i)=0.
2711	endif
2712	else
2713	proba_notrig(i)=1.
2714	random_notrig(i)=0.
2715	ale_bl_trig(i)=0.
2716	endif
2717	enddo
2718	c
2719	IF (prt_level .GE. 10) THEN
2720	print *,'proba_notrig, ale_bl_trig ',
2721	$ proba_notrig, ale_bl_trig
2722	ENDIF
2723
2724	endif !(iflag_trig_bl)
2725
2726	c-----------Statistical closure-----------
2727	if (iflag_clos_bl.ge.1) then
2728
2729	do i=1,klon
2730	alp_bl(i)=alp_bl_stat(i)
2731	enddo
2732
2733	else
2734
2735	alp_bl_stat(:)=0.
2736
2737	endif !(iflag_clos_bl)
2738
2739	IF (prt_level .GE. 10) THEN
2740	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2741	ENDIF
2742
2743	ccc fin nrlmd le 10/04/2012
2744
2745	! ----------------------------------------------------------------------
2746	! Transport de la TKE par les panaches thermiques.
2747	! FH : 2010/02/01
2748	if (iflag_pbl.eq.10) then
2749	call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2750	s rg,paprs,pbl_tke)
2751	endif
2752	! ----------------------------------------------------------------------
2753	!IM/FH: 2011/02/23
2754	! Couplage Thermiques/Emanuel seulement si T<0
2755	if (iflag_coupl==2) then
2756	print*,'Couplage Thermiques/Emanuel seulement si T<0'
2757	do i=1,klon
2758	if (t_seri(i,lmax_th(i))>273.) then
2759	Ale_bl(i)=0.
2760	endif
2761	enddo
2762	endif
2763
2764	do i=1,klon
2765	zmax_th(i)=pphi(i,lmax_th(i))/rg
2766	enddo
2767
2768	endif
2769
2770
2771	c Ajustement sec
2772	c ==============
2773
2774	! Dans le cas oÃ¹ on active les thermiques, on fait partir l'ajustement
2775	! a partir du sommet des thermiques.
2776	! Dans le cas contraire, on demarre au niveau 1.
2777
2778	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2779
2780	if(iflag_thermals.eq.0) then
2781	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2782	limbas(:)=1
2783	else
2784	limbas(:)=lmax_th(:)
2785	endif
2786
2787	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2788	! pour des test de convergence numerique.
2789	! Le nouveau ajsec est a priori mieux, meme pour le cas
2790	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2791	! non nulles numeriquement pour des mailles non concernees.
2792
2793	if (iflag_thermals.eq.0) then
2794	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2795	s , d_t_ajsb, d_q_ajsb)
2796	else
2797	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2798	s , d_t_ajsb, d_q_ajsb)
2799	endif
2800
2801	!-----------------------------------------------------------------------------------------
2802	! ajout des tendances de l'ajustement sec ou des thermiques
2803	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2804	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2805	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2806
2807	!-----------------------------------------------------------------------------------------
2808
2809	endif
2810
2811	endif
2812	c
2813	c===================================================================
2814	cIM
2815	IF (ip_ebil_phy.ge.2) THEN
2816	ztit='after dry_adjust'
2817	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2818	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2819	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2820	call diagphy(airephy,ztit,ip_ebil_phy
2821	e , zero_v, zero_v, zero_v, zero_v, zero_v
2822	e , zero_v, zero_v, zero_v, ztsol
2823	e , d_h_vcol, d_qt, d_ec
2824	s , fs_bound, fq_bound )
2825	END IF
2826
2827
2828	c-------------------------------------------------------------------------
2829	c Caclul des ratqs
2830	c-------------------------------------------------------------------------
2831
2832	c print*,'calcul des ratqs'
2833	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2834	c ----------------
2835	c on ecrase le tableau ratqsc calcule par clouds_gno
2836	if (iflag_cldcon.eq.1) then
2837	do k=1,klev
2838	do i=1,klon
2839	if(ptconv(i,k)) then
2840	ratqsc(i,k)=ratqsbas
2841	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2842	else
2843	ratqsc(i,k)=0.
2844	endif
2845	enddo
2846	enddo
2847
2848	c-----------------------------------------------------------------------
2849	c par nversion de la fonction log normale
2850	c-----------------------------------------------------------------------
2851	else if (iflag_cldcon.eq.4) then
2852	ptconvth(:,:)=.false.
2853	ratqsc(:,:)=0.
2854	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2855	call clouds_gno
2856	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2857	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2858
2859	endif
2860
2861	c ratqs stables
2862	c -------------
2863
2864	if (iflag_ratqs.eq.0) then
2865
2866	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2867	do k=1,klev
2868	do i=1, klon
2869	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2870	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2871	enddo
2872	enddo
2873
2874	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2875	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2876	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2877	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2878	! Il s'agit de differents tests dans la phase de reglage du modele
2879	! avec thermiques.
2880
2881	else if (iflag_ratqs.eq.1) then
2882
2883	do k=1,klev
2884	do i=1, klon
2885	if (pplay(i,k).ge.60000.) then
2886	ratqss(i,k)=ratqsbas
2887	else if ((pplay(i,k).ge.30000.).and.
2888	s (pplay(i,k).lt.60000.)) then
2889	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2890	s (60000.-pplay(i,k))/(60000.-30000.)
2891	else
2892	ratqss(i,k)=ratqshaut
2893	endif
2894	enddo
2895	enddo
2896
2897	else if (iflag_ratqs.eq.2) then
2898
2899	do k=1,klev
2900	do i=1, klon
2901	if (pplay(i,k).ge.60000.) then
2902	ratqss(i,k)=ratqsbas
2903	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2904	else if ((pplay(i,k).ge.30000.).and.
2905	s (pplay(i,k).lt.60000.)) then
2906	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2907	s (60000.-pplay(i,k))/(60000.-30000.)
2908	else
2909	ratqss(i,k)=ratqshaut
2910	endif
2911	enddo
2912	enddo
2913
2914	else if (iflag_ratqs==3) then
2915	do k=1,klev
2916	ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)
2917	s min( ((paprs(:,1)-pplay(:,k))/70000.)*2 , 1. )
2918	enddo
2919
2920	else if (iflag_ratqs==4) then
2921	do k=1,klev
2922	ratqss(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas)
2923	s *( tanh( (50000.-pplay(:,k))/20000.) + 1.)
2924	enddo
2925
2926	endif
2927
2928
2929
2930
2931	c ratqs final
2932	c -----------
2933
2934	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2935	s .or.iflag_cldcon.eq.4) then
2936
2937	! On ajoute une constante au ratqsc*2 pour tenir compte de
2938	! fluctuations turbulentes de petite echelle
2939
2940	do k=1,klev
2941	do i=1,klon
2942	if ((fm_therm(i,k).gt.1.e-10)) then
2943	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2944	endif
2945	enddo
2946	enddo
2947
2948	! les ratqs sont une combinaison de ratqss et ratqsc
2949	if(prt_level.ge.9)
2950	$ write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2951
2952	if (tau_ratqs>1.e-10) then
2953	facteur=exp(-pdtphys/tau_ratqs)
2954	else
2955	facteur=0.
2956	endif
2957	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2958	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2959	! FH 22/09/2009
2960	! La ligne ci-dessous faisait osciller le modele et donnait une solution
2961	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
2962	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2963	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2964	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2965	else if (iflag_cldcon<=6) then
2966	! on ne prend que le ratqs stable pour fisrtilp
2967	ratqs(:,:)=ratqss(:,:)
2968	else
2969	zfratqs1=exp(-pdtphys/10800.)
2970	zfratqs2=exp(-pdtphys/10800.)
2971	! print*,'RAPPEL RATQS 1 ',zfratqs1,zfratqs2
2972	! s ,ratqss(1,14),ratqs(1,14),ratqsc(1,14)
2973	do k=1,klev
2974	do i=1,klon
2975	if (ratqsc(i,k).gt.1.e-10) then
2976	ratqs(i,k)=ratqs(i,k)*zfratqs2
2977	s +(iflag_cldcon/100.)ratqsc(i,k)(1.-zfratqs2)
2978	endif
2979	ratqs(i,k)=min(ratqs(i,k)*zfratqs1
2980	s +ratqss(i,k)*(1.-zfratqs1),0.5)
2981	enddo
2982	enddo
2983	endif
2984
2985
2986	c
2987	c Appeler le processus de condensation a grande echelle
2988	c et le processus de precipitation
2989	c-------------------------------------------------------------------------
2990	IF (prt_level .GE.10) THEN
2991	print *,' ->fisrtilp '
2992	ENDIF
2993	c-------------------------------------------------------------------------
2994	CALL fisrtilp(dtime,paprs,pplay,
2995	. t_seri, q_seri,ptconv,ratqs,
2996	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2997	. rain_lsc, snow_lsc,
2998	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2999	. frac_impa, frac_nucl,
3000	. prfl, psfl, rhcl,
3001	. zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
3002
3003	WHERE (rain_lsc < 0) rain_lsc = 0.
3004	WHERE (snow_lsc < 0) snow_lsc = 0.
3005	!-----------------------------------------------------------------------------------------
3006	! ajout des tendances de la diffusion turbulente
3007	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
3008	!-----------------------------------------------------------------------------------------
3009	DO k = 1, klev
3010	DO i = 1, klon
3011	cldfra(i,k) = rneb(i,k)
3012	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3013	ENDDO
3014	ENDDO
3015	IF (check) THEN
3016	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
3017	WRITE(lunout,*)"apresilp=", za
3018	zx_t = 0.0
3019	za = 0.0
3020	DO i = 1, klon
3021	za = za + airephy(i)/REAL(klon)
3022	zx_t = zx_t + (rain_lsc(i)
3023	. + snow_lsc(i))*airephy(i)/REAL(klon)
3024	ENDDO
3025	zx_t = zx_t/za*dtime
3026	WRITE(lunout,*)"Precip=", zx_t
3027	ENDIF
3028	cIM
3029	IF (ip_ebil_phy.ge.2) THEN
3030	ztit='after fisrt'
3031	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3032	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3033	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3034	call diagphy(airephy,ztit,ip_ebil_phy
3035	e , zero_v, zero_v, zero_v, zero_v, zero_v
3036	e , zero_v, rain_lsc, snow_lsc, ztsol
3037	e , d_h_vcol, d_qt, d_ec
3038	s , fs_bound, fq_bound )
3039	END IF
3040
3041	if (mydebug) then
3042	call writefield_phy('u_seri',u_seri,llm)
3043	call writefield_phy('v_seri',v_seri,llm)
3044	call writefield_phy('t_seri',t_seri,llm)
3045	call writefield_phy('q_seri',q_seri,llm)
3046	endif
3047
3048	c
3049	c-------------------------------------------------------------------
3050	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3051	c-------------------------------------------------------------------
3052
3053	c 1. NUAGES CONVECTIFS
3054	c
3055	cIM cf FH
3056	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
3057	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
3058	snow_tiedtke=0.
3059	c print*,'avant calcul de la pseudo precip '
3060	c print*,'iflag_cldcon',iflag_cldcon
3061	if (iflag_cldcon.eq.-1) then
3062	rain_tiedtke=rain_con
3063	else
3064	c print*,'calcul de la pseudo precip '
3065	rain_tiedtke=0.
3066	c print*,'calcul de la pseudo precip 0'
3067	do k=1,klev
3068	do i=1,klon
3069	if (d_q_con(i,k).lt.0.) then
3070	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
3071	s *(paprs(i,k)-paprs(i,k+1))/rg
3072	endif
3073	enddo
3074	enddo
3075	endif
3076	c
3077	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3078	c
3079
3080	c Nuages diagnostiques pour Tiedtke
3081	CALL diagcld1(paprs,pplay,
3082	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
3083	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
3084	. diafra,dialiq)
3085	DO k = 1, klev
3086	DO i = 1, klon
3087	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3088	cldliq(i,k) = dialiq(i,k)
3089	cldfra(i,k) = diafra(i,k)
3090	ENDIF
3091	ENDDO
3092	ENDDO
3093
3094	ELSE IF (iflag_cldcon.ge.3) THEN
3095	c On prend pour les nuages convectifs le max du calcul de la
3096	c convection et du calcul du pas de temps precedent diminue d'un facteur
3097	c facttemps
3098	facteur = pdtphys *facttemps
3099	do k=1,klev
3100	do i=1,klon
3101	rnebcon(i,k)=rnebcon(i,k)*facteur
3102	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
3103	s then
3104	rnebcon(i,k)=rnebcon0(i,k)
3105	clwcon(i,k)=clwcon0(i,k)
3106	endif
3107	enddo
3108	enddo
3109
3110	c
3111	cjq - introduce the aerosol direct and first indirect radiative forcings
3112	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3113	IF (ok_ade.OR.ok_aie) THEN
3114	IF (.NOT. aerosol_couple)
3115	& CALL readaerosol_optic(
3116	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
3117	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
3118	& mass_solu_aero, mass_solu_aero_pi,
3119	& tau_aero, piz_aero, cg_aero,
3120	& tausum_aero, tau3d_aero)
3121	ELSE
3122	tausum_aero(:,:,:) = 0.
3123	tau_aero(:,:,:,:) = 0.
3124	piz_aero(:,:,:,:) = 0.
3125	cg_aero(:,:,:,:) = 0.
3126	ENDIF
3127
3128	cIM calcul nuages par le simulateur ISCCP
3129	c
3130	#ifdef histISCCP
3131	IF (ok_isccp) THEN
3132	c
3133	cIM lecture invtau, tautab des fichiers formattes
3134	c
3135	IF (debut) THEN
3136	c$OMP MASTER
3137	c
3138	open(99,file='tautab.formatted', FORM='FORMATTED')
3139	read(99,'(f30.20)') tautab_omp
3140	close(99)
3141	c
3142	open(99,file='invtau.formatted',form='FORMATTED')
3143	read(99,'(i10)') invtau_omp
3144
3145	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
3146	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
3147
3148	close(99)
3149	c$OMP END MASTER
3150	c$OMP BARRIER
3151	tautab=tautab_omp
3152	invtau=invtau_omp
3153	c
3154	ENDIF !debut
3155	c
3156	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
3157	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
3158	#include "calcul_simulISCCP.h"
3159	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
3160	ENDIF !ok_isccp
3161	#endif
3162
3163	c On prend la somme des fractions nuageuses et des contenus en eau
3164
3165	if (iflag_cldcon>=5) then
3166
3167	do k=1,klev
3168	ptconvth(:,k)=fm_therm(:,k+1)>0.
3169	enddo
3170
3171	if (iflag_coupl==4) then
3172
3173	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3174	! convectives et lsc dans la partie des thermiques
3175	! Le controle par iflag_coupl est peut etre provisoire.
3176	do k=1,klev
3177	do i=1,klon
3178	if (ptconv(i,k).and.ptconvth(i,k)) then
3179	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3180	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3181	else if (ptconv(i,k)) then
3182	cldfra(i,k)=rnebcon(i,k)
3183	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3184	endif
3185	enddo
3186	enddo
3187
3188	else if (iflag_coupl==5) then
3189	do k=1,klev
3190	do i=1,klon
3191	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3192	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3193	enddo
3194	enddo
3195
3196	else
3197
3198	! Si on est sur un point touche par la convection profonde et pas
3199	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3200	! de la convection profonde.
3201
3202	!IM/FH: 2011/02/23
3203	! definition des points sur lesquels ls thermiques sont actifs
3204
3205	do k=1,klev
3206	do i=1,klon
3207	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3208	cldfra(i,k)=rnebcon(i,k)
3209	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3210	endif
3211	enddo
3212	enddo
3213
3214	endif
3215
3216	else
3217
3218	! Ancienne version
3219	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3220	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3221	endif
3222
3223	ENDIF
3224
3225	! plulsc(:)=0.
3226	! do k=1,klev,-1
3227	! do i=1,klon
3228	! zzz=prfl(:,k)+psfl(:,k)
3229	! if (.not.ptconvth.zzz.gt.0.)
3230	! enddo prfl, psfl,
3231	! enddo
3232	c
3233	c 2. NUAGES STARTIFORMES
3234	c
3235	IF (ok_stratus) THEN
3236	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3237	DO k = 1, klev
3238	DO i = 1, klon
3239	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3240	cldliq(i,k) = dialiq(i,k)
3241	cldfra(i,k) = diafra(i,k)
3242	ENDIF
3243	ENDDO
3244	ENDDO
3245	ENDIF
3246	c
3247	c Precipitation totale
3248	c
3249	DO i = 1, klon
3250	rain_fall(i) = rain_con(i) + rain_lsc(i)
3251	snow_fall(i) = snow_con(i) + snow_lsc(i)
3252	ENDDO
3253	cIM
3254	IF (ip_ebil_phy.ge.2) THEN
3255	ztit="after diagcld"
3256	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3257	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3258	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3259	call diagphy(airephy,ztit,ip_ebil_phy
3260	e , zero_v, zero_v, zero_v, zero_v, zero_v
3261	e , zero_v, zero_v, zero_v, ztsol
3262	e , d_h_vcol, d_qt, d_ec
3263	s , fs_bound, fq_bound )
3264	END IF
3265	c
3266	c Calculer l'humidite relative pour diagnostique
3267	c
3268	DO k = 1, klev
3269	DO i = 1, klon
3270	zx_t = t_seri(i,k)
3271	IF (thermcep) THEN
3272	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3273	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3274	zx_qs = MIN(0.5,zx_qs)
3275	zcor = 1./(1.-retv*zx_qs)
3276	zx_qs = zx_qs*zcor
3277	ELSE
3278	IF (zx_t.LT.t_coup) THEN
3279	zx_qs = qsats(zx_t)/pplay(i,k)
3280	ELSE
3281	zx_qs = qsatl(zx_t)/pplay(i,k)
3282	ENDIF
3283	ENDIF
3284	zx_rh(i,k) = q_seri(i,k)/zx_qs
3285	zqsat(i,k)=zx_qs
3286	ENDDO
3287	ENDDO
3288
3289	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3290	c equivalente a 2m (tpote) pour diagnostique
3291	c
3292	DO i = 1, klon
3293	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3294	IF (thermcep) THEN
3295	IF(zt2m(i).LT.RTT) then
3296	Lheat=RLSTT
3297	ELSE
3298	Lheat=RLVTT
3299	ENDIF
3300	ELSE
3301	IF (zt2m(i).LT.RTT) THEN
3302	Lheat=RLSTT
3303	ELSE
3304	Lheat=RLVTT
3305	ENDIF
3306	ENDIF
3307	tpote(i) = tpot(i)*
3308	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3309	ENDDO
3310
3311	IF (type_trac == 'inca') THEN
3312	#ifdef INCA
3313	CALL VTe(VTphysiq)
3314	CALL VTb(VTinca)
3315	calday = REAL(days_elapsed + 1) + jH_cur
3316
3317	call chemtime(itap+itau_phy-1, date0, dtime)
3318	IF (config_inca == 'aero') THEN
3319	CALL AEROSOL_METEO_CALC(
3320	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
3321	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3322	END IF
3323
3324	zxsnow_dummy(:) = 0.0
3325
3326	CALL chemhook_begin (calday,
3327	$ days_elapsed+1,
3328	$ jH_cur,
3329	$ pctsrf(1,1),
3330	$ rlat,
3331	$ rlon,
3332	$ airephy,
3333	$ paprs,
3334	$ pplay,
3335	$ coefh,
3336	$ pphi,
3337	$ t_seri,
3338	$ u,
3339	$ v,
3340	$ wo(:, :, 1),
3341	$ q_seri,
3342	$ zxtsol,
3343	$ zxsnow_dummy,
3344	$ solsw,
3345	$ albsol1,
3346	$ rain_fall,
3347	$ snow_fall,
3348	$ itop_con,
3349	$ ibas_con,
3350	$ cldfra,
3351	$ iim,
3352	$ jjm,
3353	$ tr_seri,
3354	$ ftsol,
3355	$ paprs,
3356	$ cdragh,
3357	$ cdragm,
3358	$ pctsrf,
3359	$ pdtphys,
3360	$ itap)
3361
3362	CALL VTe(VTinca)
3363	CALL VTb(VTphysiq)
3364	#endif
3365	END IF !type_trac = inca
3366	c
3367	c Calculer les parametres optiques des nuages et quelques
3368	c parametres pour diagnostiques:
3369	c
3370
3371	IF (aerosol_couple) THEN
3372	mass_solu_aero(:,:) = ccm(:,:,1)
3373	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3374	END IF
3375
3376	if (ok_newmicro) then
3377	CALL newmicro (paprs, pplay,ok_newmicro,
3378	. t_seri, cldliq, cldfra, cldtau, cldemi,
3379	. cldh, cldl, cldm, cldt, cldq,
3380	. flwp, fiwp, flwc, fiwc,
3381	e ok_aie,
3382	e mass_solu_aero, mass_solu_aero_pi,
3383	e bl95_b0, bl95_b1,
3384	s cldtaupi, re, fl, ref_liq, ref_ice)
3385	else
3386	CALL nuage (paprs, pplay,
3387	. t_seri, cldliq, cldfra, cldtau, cldemi,
3388	. cldh, cldl, cldm, cldt, cldq,
3389	e ok_aie,
3390	e mass_solu_aero, mass_solu_aero_pi,
3391	e bl95_b0, bl95_b1,
3392	s cldtaupi, re, fl)
3393
3394	endif
3395	c
3396	cIM betaCRF
3397	c
3398	cldtaurad = cldtau
3399	cldemirad = cldemi
3400	c
3401	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
3402	$lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3403	c
3404	c global
3405	c
3406	DO k=1, klev
3407	DO i=1, klon
3408	if (pplay(i,k).GE.pfree) THEN
3409	beta(i,k) = beta_pbl
3410	else
3411	beta(i,k) = beta_free
3412	endif
3413	if (mskocean_beta) THEN
3414	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3415	endif
3416	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3417	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3418	ENDDO
3419	ENDDO
3420	c
3421	else
3422	c
3423	c regional
3424	c
3425	DO k=1, klev
3426	DO i=1,klon
3427	c
3428	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND.
3429	$ rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3430	if (pplay(i,k).GE.pfree) THEN
3431	beta(i,k) = beta_pbl
3432	else
3433	beta(i,k) = beta_free
3434	endif
3435	if (mskocean_beta) THEN
3436	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3437	endif
3438	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3439	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3440	endif
3441	c
3442	ENDDO
3443	ENDDO
3444	c
3445	endif
3446	c
3447	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3448	c
3449	IF (MOD(itaprad,radpas).EQ.0) THEN
3450
3451	DO i = 1, klon
3452	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
3453	. + falb1(i,is_lic) * pctsrf(i,is_lic)
3454	. + falb1(i,is_ter) * pctsrf(i,is_ter)
3455	. + falb1(i,is_sic) * pctsrf(i,is_sic)
3456	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
3457	. + falb2(i,is_lic) * pctsrf(i,is_lic)
3458	. + falb2(i,is_ter) * pctsrf(i,is_ter)
3459	. + falb2(i,is_sic) * pctsrf(i,is_sic)
3460	ENDDO
3461
3462	if (mydebug) then
3463	call writefield_phy('u_seri',u_seri,llm)
3464	call writefield_phy('v_seri',v_seri,llm)
3465	call writefield_phy('t_seri',t_seri,llm)
3466	call writefield_phy('q_seri',q_seri,llm)
3467	endif
3468
3469	IF (aerosol_couple) THEN
3470	#ifdef INCA
3471	CALL radlwsw_inca
3472	e (kdlon,kflev,dist, rmu0, fract, solaire,
3473	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
3474	e wo(:, :, 1),
3475	e cldfra, cldemirad, cldtaurad,
3476	s heat,heat0,cool,cool0,radsol,albpla,
3477	s topsw,toplw,solsw,sollw,
3478	s sollwdown,
3479	s topsw0,toplw0,solsw0,sollw0,
3480	s lwdn0, lwdn, lwup0, lwup,
3481	s swdn0, swdn, swup0, swup,
3482	e ok_ade, ok_aie,
3483	e tau_aero, piz_aero, cg_aero,
3484	s topswad_aero, solswad_aero,
3485	s topswad0_aero, solswad0_aero,
3486	s topsw_aero, topsw0_aero,
3487	s solsw_aero, solsw0_aero,
3488	e cldtaupi,
3489	s topswai_aero, solswai_aero)
3490
3491	#endif
3492	ELSE
3493	c
3494	cIM calcul radiatif pour le cas actuel
3495	c
3496	RCO2 = RCO2_act
3497	RCH4 = RCH4_act
3498	RN2O = RN2O_act
3499	RCFC11 = RCFC11_act
3500	RCFC12 = RCFC12_act
3501	c
3502	IF (prt_level .GE.10) THEN
3503	print *,' ->radlwsw, number 1 '
3504	ENDIF
3505	c
3506	CALL radlwsw
3507	e (dist, rmu0, fract,
3508	e paprs, pplay,zxtsol,albsol1, albsol2,
3509	e t_seri,q_seri,wo,
3510	e cldfra, cldemirad, cldtaurad,
3511	e ok_ade, ok_aie,
3512	e tau_aero, piz_aero, cg_aero,
3513	e cldtaupi,new_aod,
3514	e zqsat, flwc, fiwc,
3515	s heat,heat0,cool,cool0,radsol,albpla,
3516	s topsw,toplw,solsw,sollw,
3517	s sollwdown,
3518	s topsw0,toplw0,solsw0,sollw0,
3519	s lwdn0, lwdn, lwup0, lwup,
3520	s swdn0, swdn, swup0, swup,
3521	s topswad_aero, solswad_aero,
3522	s topswai_aero, solswai_aero,
3523	o topswad0_aero, solswad0_aero,
3524	o topsw_aero, topsw0_aero,
3525	o solsw_aero, solsw0_aero,
3526	o topswcf_aero, solswcf_aero)
3527
3528	c
3529	cIM 2eme calcul radiatif pour le cas perturbe ou au moins un
3530	cIM des taux doit etre different du taux actuel
3531	cIM Par defaut on a les taux perturbes egaux aux taux actuels
3532	c
3533	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.
3534	$RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR.
3535	$RCFC12_per.NE.RCFC12_act) THEN
3536	c
3537	RCO2 = RCO2_per
3538	RCH4 = RCH4_per
3539	RN2O = RN2O_per
3540	RCFC11 = RCFC11_per
3541	RCFC12 = RCFC12_per
3542	c
3543	IF (prt_level .GE.10) THEN
3544	print *,' ->radlwsw, number 2 '
3545	ENDIF
3546	c
3547	CALL radlwsw
3548	e (dist, rmu0, fract,
3549	e paprs, pplay,zxtsol,albsol1, albsol2,
3550	e t_seri,q_seri,wo,
3551	e cldfra, cldemi, cldtau,
3552	e ok_ade, ok_aie,
3553	e tau_aero, piz_aero, cg_aero,
3554	e cldtaupi,new_aod,
3555	e zqsat, flwc, fiwc,
3556	s heatp,heat0p,coolp,cool0p,radsolp,albplap,
3557	s topswp,toplwp,solswp,sollwp,
3558	s sollwdownp,
3559	s topsw0p,toplw0p,solsw0p,sollw0p,
3560	s lwdn0p, lwdnp, lwup0p, lwupp,
3561	s swdn0p, swdnp, swup0p, swupp,
3562	s topswad_aerop, solswad_aerop,
3563	s topswai_aerop, solswai_aerop,
3564	o topswad0_aerop, solswad0_aerop,
3565	o topsw_aerop, topsw0_aerop,
3566	o solsw_aerop, solsw0_aerop,
3567	o topswcf_aerop, solswcf_aerop)
3568	endif
3569	c
3570	ENDIF ! aerosol_couple
3571	itaprad = 0
3572	ENDIF ! MOD(itaprad,radpas)
3573	itaprad = itaprad + 1
3574
3575	IF (iflag_radia.eq.0) THEN
3576	IF (prt_level.ge.9) THEN
3577	PRINT *,'--------------------------------------------------'
3578	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3579	PRINT *,'>>>> heat et cool mis a zero '
3580	PRINT *,'--------------------------------------------------'
3581	END IF
3582	heat=0.
3583	cool=0.
3584	sollw=0. ! MPL 01032011
3585	solsw=0.
3586	radsol=0.
3587	END IF
3588
3589	c
3590	c Ajouter la tendance des rayonnements (tous les pas)
3591	c
3592	DO k = 1, klev
3593	DO i = 1, klon
3594	t_seri(i,k) = t_seri(i,k)
3595	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
3596	ENDDO
3597	ENDDO
3598	c
3599	if (mydebug) then
3600	call writefield_phy('u_seri',u_seri,llm)
3601	call writefield_phy('v_seri',v_seri,llm)
3602	call writefield_phy('t_seri',t_seri,llm)
3603	call writefield_phy('q_seri',q_seri,llm)
3604	endif
3605
3606	cIM
3607	IF (ip_ebil_phy.ge.2) THEN
3608	ztit='after rad'
3609	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3610	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3611	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3612	call diagphy(airephy,ztit,ip_ebil_phy
3613	e , topsw, toplw, solsw, sollw, zero_v
3614	e , zero_v, zero_v, zero_v, ztsol
3615	e , d_h_vcol, d_qt, d_ec
3616	s , fs_bound, fq_bound )
3617	END IF
3618	c
3619	c
3620	c Calculer l'hydrologie de la surface
3621	c
3622	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3623	c . agesno, ftsol,fqsurf,fsnow, ruis)
3624	c
3625
3626	c
3627	c Calculer le bilan du sol et la derive de temperature (couplage)
3628	c
3629	DO i = 1, klon
3630	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3631	c a la demande de JLD
3632	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3633	ENDDO
3634	c
3635	cmoddeblott(jan95)
3636	c Appeler le programme de parametrisation de l'orographie
3637	c a l'echelle sous-maille:
3638	c
3639	IF (prt_level .GE.10) THEN
3640	print *,' call orography ? ', ok_orodr
3641	ENDIF
3642	c
3643	IF (ok_orodr) THEN
3644	c
3645	c selection des points pour lesquels le shema est actif:
3646	igwd=0
3647	DO i=1,klon
3648	itest(i)=0
3649	c IF ((zstd(i).gt.10.0)) THEN
3650	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3651	itest(i)=1
3652	igwd=igwd+1
3653	idx(igwd)=i
3654	ENDIF
3655	ENDDO
3656	c igwdim=MAX(1,igwd)
3657	c
3658	IF (ok_strato) THEN
3659
3660	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3661	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3662	e igwd,idx,itest,
3663	e t_seri, u_seri, v_seri,
3664	s zulow, zvlow, zustrdr, zvstrdr,
3665	s d_t_oro, d_u_oro, d_v_oro)
3666
3667	ELSE
3668	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3669	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3670	e igwd,idx,itest,
3671	e t_seri, u_seri, v_seri,
3672	s zulow, zvlow, zustrdr, zvstrdr,
3673	s d_t_oro, d_u_oro, d_v_oro)
3674	ENDIF
3675	c
3676	c ajout des tendances
3677	!-----------------------------------------------------------------------------------------
3678	! ajout des tendances de la trainee de l'orographie
3679	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3680	!-----------------------------------------------------------------------------------------
3681	c
3682	ENDIF ! fin de test sur ok_orodr
3683	c
3684	if (mydebug) then
3685	call writefield_phy('u_seri',u_seri,llm)
3686	call writefield_phy('v_seri',v_seri,llm)
3687	call writefield_phy('t_seri',t_seri,llm)
3688	call writefield_phy('q_seri',q_seri,llm)
3689	endif
3690
3691	IF (ok_orolf) THEN
3692	c
3693	c selection des points pour lesquels le shema est actif:
3694	igwd=0
3695	DO i=1,klon
3696	itest(i)=0
3697	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3698	itest(i)=1
3699	igwd=igwd+1
3700	idx(igwd)=i
3701	ENDIF
3702	ENDDO
3703	c igwdim=MAX(1,igwd)
3704	c
3705	IF (ok_strato) THEN
3706
3707	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3708	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3709	e igwd,idx,itest,
3710	e t_seri, u_seri, v_seri,
3711	s zulow, zvlow, zustrli, zvstrli,
3712	s d_t_lif, d_u_lif, d_v_lif )
3713
3714	ELSE
3715	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3716	e rlat,zmea,zstd,zpic,
3717	e itest,
3718	e t_seri, u_seri, v_seri,
3719	s zulow, zvlow, zustrli, zvstrli,
3720	s d_t_lif, d_u_lif, d_v_lif)
3721	ENDIF
3722	c
3723	!-----------------------------------------------------------------------------------------
3724	! ajout des tendances de la portance de l'orographie
3725	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3726	!-----------------------------------------------------------------------------------------
3727	c
3728	ENDIF ! fin de test sur ok_orolf
3729	C HINES GWD PARAMETRIZATION
3730
3731	IF (ok_hines) then
3732
3733	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3734	i rlat,t_seri,u_seri,v_seri,
3735	o zustrhi,zvstrhi,
3736	o d_t_hin, d_u_hin, d_v_hin)
3737	c
3738	c ajout des tendances
3739	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
3740
3741	ENDIF
3742	c
3743
3744	c
3745	cIM cf. FLott BEG
3746	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3747
3748	if (mydebug) then
3749	call writefield_phy('u_seri',u_seri,llm)
3750	call writefield_phy('v_seri',v_seri,llm)
3751	call writefield_phy('t_seri',t_seri,llm)
3752	call writefield_phy('q_seri',q_seri,llm)
3753	endif
3754
3755	DO i = 1, klon
3756	zustrph(i)=0.
3757	zvstrph(i)=0.
3758	ENDDO
3759	DO k = 1, klev
3760	DO i = 1, klon
3761	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3762	c (paprs(i,k)-paprs(i,k+1))/rg
3763	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3764	c (paprs(i,k)-paprs(i,k+1))/rg
3765	ENDDO
3766	ENDDO
3767	c
3768	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3769	c
3770	IF (is_sequential .and. ok_orodr) THEN
3771	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3772	C ra,rg,romega,
3773	C rlat,rlon,pphis,
3774	C zustrdr,zustrli,zustrph,
3775	C zvstrdr,zvstrli,zvstrph,
3776	C paprs,u,v,
3777	C aam, torsfc)
3778	ENDIF
3779	cIM cf. FLott END
3780	cIM
3781	IF (ip_ebil_phy.ge.2) THEN
3782	ztit='after orography'
3783	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3784	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3785	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3786	call diagphy(airephy,ztit,ip_ebil_phy
3787	e , zero_v, zero_v, zero_v, zero_v, zero_v
3788	e , zero_v, zero_v, zero_v, ztsol
3789	e , d_h_vcol, d_qt, d_ec
3790	s , fs_bound, fq_bound )
3791	END IF
3792	c
3793	c
3794	!====================================================================
3795	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3796	!====================================================================
3797	! Abderrahmane 24.08.09
3798
3799	IF (ok_cosp) THEN
3800	! adeclarer
3801	#ifdef CPP_COSP
3802	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3803
3804	print*,'freq_cosp',freq_cosp
3805	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3806	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3807	! s ref_liq,ref_ice
3808	call phys_cosp(itap,dtime,freq_cosp,
3809	$ ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
3810	$ ecrit_mth,ecrit_day,ecrit_hf,
3811	$ klon,klev,rlon,rlat,presnivs,overlap,
3812	$ ref_liq,ref_ice,
3813	$ pctsrf(:,is_ter)+pctsrf(:,is_lic),
3814	$ zu10m,zv10m,pphis,
3815	$ zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3816	$ qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3817	$ prfl(:,1:klev),psfl(:,1:klev),
3818	$ pmflxr(:,1:klev),pmflxs(:,1:klev),
3819	$ mr_ozone,cldtaurad, cldemirad)
3820
3821	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3822	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3823	! M clMISR,
3824	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3825	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3826
3827	ENDIF
3828
3829	#endif
3830	ENDIF !ok_cosp
3831	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3832	cAA
3833	cAA Installation de l'interface online-offline pour traceurs
3834	cAA
3835	c====================================================================
3836	c Calcul des tendances traceurs
3837	c====================================================================
3838	C
3839
3840	IF (type_trac=='repr') THEN
3841	sh_in(:,:) = q_seri(:,:)
3842	ELSE
3843	sh_in(:,:) = qx(:,:,ivap)
3844	END IF
3845
3846	call phytrac (
3847	I itap, days_elapsed+1, jH_cur, debut,
3848	I lafin, dtime, u, v, t,
3849	I paprs, pplay, pmfu, pmfd,
3850	I pen_u, pde_u, pen_d, pde_d,
3851	I cdragh, coefh, fm_therm, entr_therm,
3852	I u1, v1, ftsol, pctsrf,
3853	I rlat, frac_impa, frac_nucl,rlon,
3854	I presnivs, pphis, pphi, albsol1,
3855	I sh_in, rhcl, cldfra, rneb,
3856	I diafra, cldliq, itop_con, ibas_con,
3857	I pmflxr, pmflxs, prfl, psfl,
3858	I da, phi, mp, upwd,
3859	I dnwd, aerosol_couple, flxmass_w,
3860	I tau_aero, piz_aero, cg_aero, ccm,
3861	I rfname,
3862	O tr_seri)
3863
3864	IF (offline) THEN
3865
3866	IF (prt_level.ge.9)
3867	$ print*,'Attention on met a 0 les thermiques pour phystoke'
3868	call phystokenc (
3869	I nlon,klev,pdtphys,rlon,rlat,
3870	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3871	I fm_therm,entr_therm,
3872	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3873	I frac_impa, frac_nucl,
3874	I pphis,airephy,dtime,itap,
3875	I qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3876
3877
3878	ENDIF
3879
3880	c
3881	c Calculer le transport de l'eau et de l'energie (diagnostique)
3882	c
3883	CALL transp (paprs,zxtsol,
3884	e t_seri, q_seri, u_seri, v_seri, zphi,
3885	s ve, vq, ue, uq)
3886	c
3887	cIM global posePB BEG
3888	IF(1.EQ.0) THEN
3889	c
3890	CALL transp_lay (paprs,zxtsol,
3891	e t_seri, q_seri, u_seri, v_seri, zphi,
3892	s ve_lay, vq_lay, ue_lay, uq_lay)
3893	c
3894	ENDIF !(1.EQ.0) THEN
3895	cIM global posePB END
3896	c Accumuler les variables a stocker dans les fichiers histoire:
3897	c
3898	c+jld ec_conser
3899	DO k = 1, klev
3900	DO i = 1, klon
3901	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3902	d_t_ec(i,k)=0.5/ZRCPD
3903	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3904	ENDDO
3905	ENDDO
3906
3907	DO k = 1, klev
3908	DO i = 1, klon
3909	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3910	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3911	END DO
3912	END DO
3913	c-jld ec_conser
3914	cIM
3915	IF (ip_ebil_phy.ge.1) THEN
3916	ztit='after physic'
3917	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3918	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3919	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3920	C Comme les tendances de la physique sont ajoute dans la dynamique,
3921	C on devrait avoir que la variation d'entalpie par la dynamique
3922	C est egale a la variation de la physique au pas de temps precedent.
3923	C Donc la somme de ces 2 variations devrait etre nulle.
3924
3925	call diagphy(airephy,ztit,ip_ebil_phy
3926	e , topsw, toplw, solsw, sollw, sens
3927	e , evap, rain_fall, snow_fall, ztsol
3928	e , d_h_vcol, d_qt, d_ec
3929	s , fs_bound, fq_bound )
3930	C
3931	d_h_vcol_phy=d_h_vcol
3932	C
3933	END IF
3934	C
3935	c=======================================================================
3936	c SORTIES
3937	c=======================================================================
3938
3939	cIM Interpolation sur les niveaux de pression du NMC
3940	c -------------------------------------------------
3941	c
3942	#include "calcul_STDlev.h"
3943	twriteSTD(:,:,1)=tsumSTD(:,:,1)
3944	qwriteSTD(:,:,1)=qsumSTD(:,:,1)
3945	rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
3946	phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
3947	uwriteSTD(:,:,1)=usumSTD(:,:,1)
3948	vwriteSTD(:,:,1)=vsumSTD(:,:,1)
3949	wwriteSTD(:,:,1)=wsumSTD(:,:,1)
3950
3951	twriteSTD(:,:,2)=tsumSTD(:,:,2)
3952	qwriteSTD(:,:,2)=qsumSTD(:,:,2)
3953	rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
3954	phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
3955	uwriteSTD(:,:,2)=usumSTD(:,:,2)
3956	vwriteSTD(:,:,2)=vsumSTD(:,:,2)
3957	wwriteSTD(:,:,2)=wsumSTD(:,:,2)
3958
3959	twriteSTD(:,:,3)=tlevSTD(:,:)
3960	qwriteSTD(:,:,3)=qlevSTD(:,:)
3961	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3962	phiwriteSTD(:,:,3)=philevSTD(:,:)
3963	uwriteSTD(:,:,3)=ulevSTD(:,:)
3964	vwriteSTD(:,:,3)=vlevSTD(:,:)
3965	wwriteSTD(:,:,3)=wlevSTD(:,:)
3966
3967	twriteSTD(:,:,4)=tlevSTD(:,:)
3968	qwriteSTD(:,:,4)=qlevSTD(:,:)
3969	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3970	phiwriteSTD(:,:,4)=philevSTD(:,:)
3971	uwriteSTD(:,:,4)=ulevSTD(:,:)
3972	vwriteSTD(:,:,4)=vlevSTD(:,:)
3973	wwriteSTD(:,:,4)=wlevSTD(:,:)
3974	c
3975	cIM initialisation 5eme fichier de sortie
3976	twriteSTD(:,:,5)=tlevSTD(:,:)
3977	qwriteSTD(:,:,5)=qlevSTD(:,:)
3978	rhwriteSTD(:,:,5)=rhlevSTD(:,:)
3979	phiwriteSTD(:,:,5)=philevSTD(:,:)
3980	uwriteSTD(:,:,5)=ulevSTD(:,:)
3981	vwriteSTD(:,:,5)=vlevSTD(:,:)
3982	wwriteSTD(:,:,5)=wlevSTD(:,:)
3983	c
3984	cIM initialisation 6eme fichier de sortie
3985	twriteSTD(:,:,6)=tlevSTD(:,:)
3986	qwriteSTD(:,:,6)=qlevSTD(:,:)
3987	rhwriteSTD(:,:,6)=rhlevSTD(:,:)
3988	phiwriteSTD(:,:,6)=philevSTD(:,:)
3989	uwriteSTD(:,:,6)=ulevSTD(:,:)
3990	vwriteSTD(:,:,6)=vlevSTD(:,:)
3991	wwriteSTD(:,:,6)=wlevSTD(:,:)
3992	cIM for NMC files
3993	DO n=1, nlevSTD3
3994	DO k=1, nlevSTD
3995	if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
3996	twriteSTD3(:,n)=tlevSTD(:,k)
3997	qwriteSTD3(:,n)=qlevSTD(:,k)
3998	rhwriteSTD3(:,n)=rhlevSTD(:,k)
3999	phiwriteSTD3(:,n)=philevSTD(:,k)
4000	uwriteSTD3(:,n)=ulevSTD(:,k)
4001	vwriteSTD3(:,n)=vlevSTD(:,k)
4002	wwriteSTD3(:,n)=wlevSTD(:,k)
4003	endif !rlevSTD3(n).EQ.rlevSTD(k)
4004	ENDDO
4005	ENDDO
4006	c
4007	DO n=1, nlevSTD8
4008	DO k=1, nlevSTD
4009	if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
4010	tnondefSTD8(:,n)=tnondef(:,k,2)
4011	twriteSTD8(:,n)=tsumSTD(:,k,2)
4012	qwriteSTD8(:,n)=qsumSTD(:,k,2)
4013	rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
4014	phiwriteSTD8(:,n)=phisumSTD(:,k,2)
4015	uwriteSTD8(:,n)=usumSTD(:,k,2)
4016	vwriteSTD8(:,n)=vsumSTD(:,k,2)
4017	wwriteSTD8(:,n)=wsumSTD(:,k,2)
4018	endif !rlevSTD8(n).EQ.rlevSTD(k)
4019	ENDDO
4020	ENDDO
4021	c
4022	c slp sea level pressure
4023	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
4024	c
4025	ccc prw = eau precipitable
4026	DO i = 1, klon
4027	prw(i) = 0.
4028	DO k = 1, klev
4029	prw(i) = prw(i) +
4030	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
4031	ENDDO
4032	ENDDO
4033	c
4034	cIM initialisation + calculs divers diag AMIP2
4035	c
4036	#include "calcul_divers.h"
4037	c
4038	IF (type_trac == 'inca') THEN
4039	#ifdef INCA
4040	CALL VTe(VTphysiq)
4041	CALL VTb(VTinca)
4042
4043	CALL chemhook_end (
4044	$ dtime,
4045	$ pplay,
4046	$ t_seri,
4047	$ tr_seri,
4048	$ nbtr,
4049	$ paprs,
4050	$ q_seri,
4051	$ airephy,
4052	$ pphi,
4053	$ pphis,
4054	$ zx_rh)
4055
4056	CALL VTe(VTinca)
4057	CALL VTb(VTphysiq)
4058	#endif
4059	END IF
4060
4061	c=============================================================
4062	c
4063	c Convertir les incrementations en tendances
4064	c
4065	IF (prt_level .GE.10) THEN
4066	print *,'Convertir les incrementations en tendances '
4067	ENDIF
4068	c
4069	if (mydebug) then
4070	call writefield_phy('u_seri',u_seri,llm)
4071	call writefield_phy('v_seri',v_seri,llm)
4072	call writefield_phy('t_seri',t_seri,llm)
4073	call writefield_phy('q_seri',q_seri,llm)
4074	endif
4075
4076	DO k = 1, klev
4077	DO i = 1, klon
4078	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4079	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4080	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4081	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4082	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4083	ENDDO
4084	ENDDO
4085	c
4086	IF (nqtot.GE.3) THEN
4087	DO iq = 3, nqtot
4088	DO k = 1, klev
4089	DO i = 1, klon
4090	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4091	ENDDO
4092	ENDDO
4093	ENDDO
4094	ENDIF
4095	c
4096	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4097	cIM global posePB#include "write_bilKP_ins.h"
4098	cIM global posePB#include "write_bilKP_ave.h"
4099	c
4100
4101	c Sauvegarder les valeurs de t et q a la fin de la physique:
4102	c
4103	DO k = 1, klev
4104	DO i = 1, klon
4105	u_ancien(i,k) = u_seri(i,k)
4106	v_ancien(i,k) = v_seri(i,k)
4107	t_ancien(i,k) = t_seri(i,k)
4108	q_ancien(i,k) = q_seri(i,k)
4109	ENDDO
4110	ENDDO
4111	c
4112	!==========================================================================
4113	! Sorties des tendances pour un point particulier
4114	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4115	! pour le debug
4116	! La valeur de igout est attribuee plus haut dans le programme
4117	!==========================================================================
4118
4119	if (prt_level.ge.1) then
4120	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4121	write(lunout,*)
4122	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4123	write(lunout,*)
4124	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
4125	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
4126	s pctsrf(igout,is_sic)
4127	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4128	do k=1,klev
4129	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
4130	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
4131	s d_t_eva(igout,k)
4132	enddo
4133	write(lunout,*) 'cool,heat'
4134	do k=1,klev
4135	write(lunout,*) cool(igout,k),heat(igout,k)
4136	enddo
4137
4138	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4139	do k=1,klev
4140	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
4141	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4142	enddo
4143
4144	write(lunout,*) 'd_ps ',d_ps(igout)
4145	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4146	do k=1,klev
4147	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
4148	s d_qx(igout,k,1),d_qx(igout,k,2)
4149	enddo
4150	endif
4151
4152	!==========================================================================
4153
4154	c============================================================
4155	c Calcul de la temperature potentielle
4156	c============================================================
4157	DO k = 1, klev
4158	DO i = 1, klon
4159	cJYG/IM theta en debut du pas de temps
4160	cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4161	cJYG/IM theta en fin de pas de temps de physique
4162	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4163	ENDDO
4164	ENDDO
4165	c
4166
4167	c 22.03.04 BEG
4168	c=============================================================
4169	c Ecriture des sorties
4170	c=============================================================
4171	#ifdef CPP_IOIPSL
4172
4173	c Recupere des varibles calcule dans differents modules
4174	c pour ecriture dans histxxx.nc
4175
4176	! Get some variables from module fonte_neige_mod
4177	CALL fonte_neige_get_vars(pctsrf,
4178	. zxfqcalving, zxfqfonte, zxffonte)
4179
4180
4181
4182	c=============================================================
4183	! Separation entre thermiques et non thermiques dans les sorties
4184	! de fisrtilp
4185	c=============================================================
4186
4187	if (iflag_thermals>=1) then
4188	d_t_lscth=0.
4189	d_t_lscst=0.
4190	d_q_lscth=0.
4191	d_q_lscst=0.
4192	do k=1,klev
4193	do i=1,klon
4194	if (ptconvth(i,k)) then
4195	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4196	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4197	else
4198	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4199	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4200	endif
4201	enddo
4202	enddo
4203
4204	do i=1,klon
4205	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4206	plul_th(i)=prfl(i,1)+psfl(i,1)
4207	enddo
4208	endif
4209
4210
4211	#include "phys_output_write.h"
4212
4213	#ifdef histISCCP
4214	#include "write_histISCCP.h"
4215	#endif
4216
4217	#ifdef histNMC
4218	#include "write_histhfNMC.h"
4219	#include "write_histdayNMC.h"
4220	#include "write_histmthNMC.h"
4221	#endif
4222
4223	#include "write_histday_seri.h"
4224
4225	#include "write_paramLMDZ_phy.h"
4226
4227	#endif
4228
4229	c 22.03.04 END
4230	c
4231	c====================================================================
4232	c Si c'est la fin, il faut conserver l'etat de redemarrage
4233	c====================================================================
4234	c
4235
4236	c -----------------------------------------------------------------
4237	c WSTATS: Saving statistics
4238	c -----------------------------------------------------------------
4239	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
4240	c which can later be used to make the statistic files of the run:
4241	c "stats") only possible in 3D runs !
4242
4243
4244	IF (callstats) THEN
4245
4246	call wstats(klon,o_psol%name,"Surface pressure","Pa"
4247	& ,2,paprs(:,1))
4248	call wstats(klon,o_tsol%name,"Surface temperature","K",
4249	& 2,zxtsol)
4250	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4251	call wstats(klon,o_precip%name,"Precip Totale liq+sol",
4252	& "kg/(s*m2)",2,zx_tmp_fi2d)
4253	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4254	call wstats(klon,o_plul%name,"Large-scale Precip",
4255	& "kg/(s*m2)",2,zx_tmp_fi2d)
4256	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4257	call wstats(klon,o_pluc%name,"Convective Precip",
4258	& "kg/(s*m2)",2,zx_tmp_fi2d)
4259	call wstats(klon,o_sols%name,"Solar rad. at surf.",
4260	& "W/m2",2,solsw)
4261	call wstats(klon,o_soll%name,"IR rad. at surf.",
4262	& "W/m2",2,sollw)
4263	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4264	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA",
4265	& "W/m2",2,zx_tmp_fi2d)
4266
4267
4268
4269	call wstats(klon,o_temp%name,"Air temperature","K",
4270	& 3,t_seri)
4271	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1",
4272	& 3,u_seri)
4273	call wstats(klon,o_vitv%name,"Meridional wind",
4274	& "m.s-1",3,v_seri)
4275	call wstats(klon,o_vitw%name,"Vertical wind",
4276	& "m.s-1",3,omega)
4277	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg",
4278	& 3,q_seri)
4279
4280
4281
4282	IF(lafin) THEN
4283	write (,) "Writing stats..."
4284	call mkstats(ierr)
4285	ENDIF
4286
4287	ENDIF !if callstats
4288
4289
4290	IF (lafin) THEN
4291	itau_phy = itau_phy + itap
4292	CALL phyredem ("restartphy.nc")
4293	! open(97,form="unformatted",file="finbin")
4294	! write(97) u_seri,v_seri,t_seri,q_seri
4295	! close(97)
4296	C$OMP MASTER
4297	if (read_climoz >= 1) then
4298	if (is_mpi_root) then
4299	call nf95_close(ncid_climoz)
4300	end if
4301	deallocate(press_climoz) ! pointer
4302	end if
4303	C$OMP END MASTER
4304	ENDIF
4305
4306	! first=.false.
4307
4308	RETURN
4309	END
4310	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4311	IMPLICIT none
4312	c
4313	c Calculer et imprimer l'eau totale. A utiliser pour verifier
4314	c la conservation de l'eau
4315	c
4316	#include "YOMCST.h"
4317	INTEGER klon,klev
4318	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4319	REAL aire(klon)
4320	REAL qtotal, zx, qcheck
4321	INTEGER i, k
4322	c
4323	zx = 0.0
4324	DO i = 1, klon
4325	zx = zx + aire(i)
4326	ENDDO
4327	qtotal = 0.0
4328	DO k = 1, klev
4329	DO i = 1, klon
4330	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
4331	. *(paprs(i,k)-paprs(i,k+1))/RG
4332	ENDDO
4333	ENDDO
4334	c
4335	qcheck = qtotal/zx
4336	c
4337	RETURN
4338	END
4339	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4340	IMPLICIT none
4341	c
4342	c Tranformer une variable de la grille physique a
4343	c la grille d'ecriture
4344	c
4345	INTEGER nfield,nlon,iim,jjmp1, jjm
4346	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4347	c
4348	INTEGER i, n, ig
4349	c
4350	jjm = jjmp1 - 1
4351	DO n = 1, nfield
4352	DO i=1,iim
4353	ecrit(i,n) = fi(1,n)
4354	ecrit(i+jjm*iim,n) = fi(nlon,n)
4355	ENDDO
4356	DO ig = 1, nlon - 2
4357	ecrit(iim+ig,n) = fi(1+ig,n)
4358	ENDDO
4359	ENDDO
4360	RETURN
4361	END

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

Download in other formats:

Original Format