Context Navigation

physiq.F @ 1813

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

transform phytrac into a module, in order to pass some variables

(tracer tendencies) to the standard physiq ouput codes.

Correct a (big) bug in the call to phytrac.
Include w*, and ALEs in the call to phytrac and traclmdz.

physiq.F

Bug correction in the call of phytrac from the physics u10m,v10m, ustar -> zu10m, zv10m, zustar

phytrac.F90 -> phytrac_mod.F90

Tranformation of routine phytrac into a module phytrac_mod, in order to tranfer the tracer tendencies from phytrac to

phys_output...

Inclusion of w*, Ale bl/wake in the call to phytrac and traclmdz

(to be used for dust emmission)

by respectively, vertical diffusion, thermal plumes and convection

desactivation of ini_histrac.h and write_histrac.h
USE phys_output_mod removed since it was creating a circular

dependency

between phytrac_mod and phys_output_mod.
So the automatic computation of ecrit_tra is desactivated

ini_histrac.h and write_histrac.h

Descactivated in phytrac but kept for backard compatibility
couchelimite -> iflag_vdf_trac>0

phys_output_ctrlout_mod.F90

New variables : o_dtr_vdf, o_dtr_the ... for output of tracer tendencies

phys_output_mod.F90

Default definition for these new output variables.

phys_output_write_F90.h

disapears, included directly in phys_output_write_mod.F90

phys_output_write_mod.F90

writing of the tracer tendencies

phys_state_var_mod.F90

New declaration (wstar)

traclmdz_mod.F90

Inclusion of w*, Ale bl/wake in the call to phytrac and traclmdz

(to be used for dust emmission)

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

File size: 133.1 KB

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

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

Download in other formats:

Original Format