Context Navigation

physiq.F @ 1885

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

Création d'une nouvelle testing:

merge des modifications du trunk entre r1796 et r1860

New testing version

merged modifications between r1796 and r1860 from the trunk

i.e.
svn merge -r1796:1860 http://svn.lmd.jussieu.fr/LMDZ/LMDZ5/trunk

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

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

Download in other formats:

Original Format