Context Navigation

physiq.F @ 1764

Last change on this file since 1764 was 1764, checked in by Laurent Fairhead, 11 years ago

Inclusion d'une routine qui lit des champs d'aérosols stratosphériques
mensuels, prescrit des propriétés optiques et modifie le rayonnement en
conséquence. Pour le moment, seule l'interaction avec le rayonnement ondes
courtes est pris en compte. Les fichiers d'input doivent être au format des
fichiers de sortie. Contrôlé par la variable logique: flag_aerosol_strat
(false par défaut dans DefLists?/config.def)

Boucher

A new routine has been added to the code that reads in monthly stratospheric
aerosols, prescribes optical properties and modifies radiation accordingly.
Presently, only the interaction with short wave radiation is taken into account.
Input files must be formatted as are the aerosol output fields. Control is by
the logical flag: flag_aerosol_strat (which is false by default and included
DefLists?/config.def)

Boucher

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

File size: 136.0 KB

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

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

Download in other formats:

Original Format