Context Navigation

physiq.F @ 1638

Last change on this file since 1638 was 1638, checked in by idelkadi, 12 years ago
Introduction du declenchement stochastique de la convection
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 142.0 KB

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

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

Download in other formats:

Original Format