Context Navigation

physiq.F @ 1712

Last change on this file since 1712 was 1712, checked in by idelkadi, 11 years ago
Corrections dans newmicro.F pour completer les modifications de O. Boucher sur ok_ade/ok_aie avec correction des diagnostiques CMIP5 : Rajout d'un nouveau flag ok_cdnc (ok cloud droplet number concentration) Dans le cas sans aérosols, nous avons flag_aerosol=0, ok_cdnc=n, ok_ade=n et ok_aie=n Dans le cas avec aérosols, nous avons flag_aerosol=6, ok_cdnc=y, ok_ade=y et ok_aie=y
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 132.5 KB

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

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

Download in other formats:

Original Format