Context Navigation

physiq.F @ 1352

Last change on this file since 1352 was 1352, checked in by musat, 14 years ago

Add 3 output files for standard pressure levels AR5 exercice and flags
to manage their computation and output frequencies
histhfNMC.nc with 3 standard pressure levels
histdayNMC.nc with 8 (or may have 17) standard pressure levels
histmthNMC.nc with 17 standard pressure levels
Add 3 flags in the physiq.def file: freq_calNMC(3), freq_outNMC(3) and lev_histdayNMC
freq_calNMC(3) : computation frequency of variables on standard pressure levels

and by default has fallowing values (in fact physics' time step dtime)

freq_calNMC(1)=900.
freq_calNMC(2)=900.
freq_calNMC(3)=900.
freq_outNMC(3) : output frequency of variables on standard pressure levels

with following default values

freq_out(1) = 2592000. (30 days)
freq_out(2) = 86400. (1 day)
freq_out(3) = 21600. (6 hours)
lev_histdayNMC is 8 by default but may be switched to 17 (if we need more levels for a particular run)
IM

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

File size: 120.0 KB

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

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

Download in other formats:

Original Format