Context Navigation

physiq.F @ 1392

Last change on this file since 1392 was 1392, checked in by musat, 14 years ago
Variables (read_climoz, ncid_climoz, press_climoz) are "OpenMP shared" and do not need $OMP THREADPRIVATE. LG/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 123.4 KB

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

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

Download in other formats:

Original Format