Context Navigation

physiq.F @ 1735

Last change on this file since 1735 was 1735, checked in by idelkadi, 11 years ago
Corrections dans le cas des simulations de sensibilite ou l'on desactive (partout ou seulement dans la couche limite) les effets radiatifs des nuages. Les parametres de controle de ces simulations sont contenus dans le fichier beta_crf.data
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 132.9 KB

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

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

Download in other formats:

Original Format