Context Navigation

physiq.F @ 1671

Last change on this file since 1671 was 1670, checked in by idelkadi, 12 years ago

Modifications for inclusion of chimere dust emission module :
u* is passed from the boundary layer parameterization to the physics
main routine (physiq.F) and then to phytrac, traclmdz and change_srf_frac.
The interface of traclmdz is enriched with 4 other variables.
Also u* and the vertically cumulated amount of tracers is added in the
outputs.

Modifications pour l'inclusion du module d'émission de poussière de Chimere :
u* est passé depuis la couche limite vers le programme principal de la
physique (physiq.F) et ensuite à phytrac, traclmdz et change_srf_frac.
L'interface de traclmdz est enrichie avec 4 autres variables.
Les variables u* et les cumuls verticaux des traceurs sont ajoutés
dans les sorties.

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

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

Download in other formats:

Original Format