Context Navigation

physiq.F @ 1753

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

Concerns energy conservation.

The source terms of the TKE prognostic equation are diagnosed
from the tendencies (du, dv, dT) associated with subrgrid scale
motions and treated as an additional heat source.
Controled by a new key, iflag_ener_conserv :

0 no conservation

-1 old adhoc correction for kinetic E only (used for CMIP5)

1 conservation
101 conversion from kinetic to heat only
110 conversion from potential to heat only

iflag_ener_conserv=-1 kept as default value for a test period
iflag_ener_conserv=1 is the advisable value
Concerns clesphys.h, and conf_phys.F90
New routine : ener_conserv.F90, called by physic.
New outputs :
bils_ec, contribution to the energy budget of the column of the

additional heat source (in W/m2)

bils_kinetic : change kinetic energy of the column in physics (W/m2)
bils_enthalp : idem for the total column enthalphy
bils_latent : idem for latent heat
Modified files : clesphys.h, conf_phys_m.F90, physiq.F,

phys_output_mod.F90, phys_output_var_mod.F90, phys_output_write.h,
ener_conserv.F90

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

File size: 135.3 KB

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

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

Download in other formats:

Original Format