Context Navigation

physiq.F @ 1565

Last change on this file since 1565 was 1565, checked in by jghattas, 13 years ago

Added interface with chemestry model REPROBUS :

Compile LMDZ together with Reprobus code (dependecies in both directions) and cpp key REPROBUS :

./makelmdz_fcm -ext_src my_path_to_reprobus -cpp REPROBUS ...

For running, add type_trac=repr in run.def.

/Marion Marchand, JG

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

File size: 136.8 KB

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

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

Download in other formats:

Original Format