Context Navigation

physiq.F @ 1456

Last change on this file since 1456 was 1456, checked in by musat, 14 years ago
phyetat0, phyredem: correction dimension verticale pbl_tke: pbl_tke(:,1:klev+1,:) physiq: pour pouvoir fixer la longitude solaire avec la nouvelle orbite JYG/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 123.0 KB

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

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

Download in other formats:

Original Format