Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

physiq.F @ 5426

Last change on this file since 5426 was 1634, checked in by Laurent Fairhead, 12 years ago

Version du code LMDZ utilisée dans la configuration IPSLCM5B
de l'exercice CMIP5/AR5

Actual version of the code used for the IPSLCM5B configuration
of the CMIP5/AR5 exercise

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

File size: 136.2 KB

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

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

Download in other formats:

Original Format