Context Navigation

physiq.F @ 1350

Last change on this file since 1350 was 1350, 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: 122.3 KB

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

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

Download in other formats:

Original Format