Context Navigation

physiq.F @ 1743

Last change on this file since 1743 was 1742, checked in by idelkadi, 11 years ago

1- Inclusion des developpements de la these de Romain Pilon sur le
lessivage des aerosols :

a/ par les pluies convectives (modifs cv30_routines et cv3_routines pour

sortir les champs nécessaires au calcul off-line ; modif cvltr)

b/ par les pluies stratiformes (modifs phytrac et introduction

lsc_scav).

2- Choix entre plusieurs schemas pour les pluies stratiformes, commande
par iflag_lscav.

3- Quelques corrections dans la convection "Nouvelle Physique" pour
assurer la conservation des traceurs (cv3p1_mixing et cva_driver) (travail
de Robin Locatelli).

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

File size: 135.2 KB

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

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

Download in other formats:

Original Format