Context Navigation

physiq.F @ 3671

Last change on this file since 3671 was 1493, checked in by Laurent Fairhead, 14 years ago

Few more fields need to be set to zero when radiation is not activated

Mise à zéro de champs supplémentaires quand le rayonnement est désactivé

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

File size: 125.8 KB

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

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

Download in other formats:

Original Format