Context Navigation

physiq.F @ 1286

Last change on this file since 1286 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 117.4 KB

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

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

Download in other formats:

Original Format