Context Navigation

physiq.F @ 1454

Last change on this file since 1454 was 1454, checked in by Laurent Fairhead, 14 years ago

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434

Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

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

File size: 123.1 KB

Line
1	! $Id: physiq.F 1454 2010-11-18 12:01:24Z fairhead $
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	cIM 180608
1253	c pmflxr=0.
1254	c pmflxs=0.
1255
1256	itau_con=0
1257	first=.false.
1258
1259	endif ! first
1260
1261	modname = 'physiq'
1262	cIM
1263	IF (ip_ebil_phy.ge.1) THEN
1264	DO i=1,klon
1265	zero_v(i)=0.
1266	END DO
1267	END IF
1268	ok_sync=.TRUE.
1269
1270	IF (debut) THEN
1271	CALL suphel ! initialiser constantes et parametres phys.
1272	ENDIF
1273
1274	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1275
1276
1277	c======================================================================
1278	! Gestion calendrier : mise a jour du module phys_cal_mod
1279	!
1280	c CALL phys_cal_update(jD_cur,jH_cur)
1281
1282	c
1283	c Si c'est le debut, il faut initialiser plusieurs choses
1284	c ********
1285	c
1286	IF (debut) THEN
1287	!rv
1288	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1289	cde la convection a partir des caracteristiques du thermique
1290	wght_th(:,:)=1.
1291	lalim_conv(:)=1
1292	cRC
1293	u10m(:,:)=0.
1294	v10m(:,:)=0.
1295	rain_con(:)=0.
1296	snow_con(:)=0.
1297	topswai(:)=0.
1298	topswad(:)=0.
1299	solswai(:)=0.
1300	solswad(:)=0.
1301
1302	lambda_th(:,:)=0.
1303	wmax_th(:)=0.
1304	tau_overturning_th(:)=0.
1305
1306	IF (config_inca /= 'none') THEN
1307	! jg : initialisation jusqu'au ces variables sont dans restart
1308	ccm(:,:,:) = 0.
1309	tau_aero(:,:,:,:) = 0.
1310	piz_aero(:,:,:,:) = 0.
1311	cg_aero(:,:,:,:) = 0.
1312	END IF
1313
1314	rnebcon0(:,:) = 0.0
1315	clwcon0(:,:) = 0.0
1316	rnebcon(:,:) = 0.0
1317	clwcon(:,:) = 0.0
1318
1319	cIM
1320	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1321	c
1322	print*,'iflag_coupl,iflag_clos,iflag_wake',
1323	. iflag_coupl,iflag_clos,iflag_wake
1324	print*,'CYCLE_DIURNE', cycle_diurne
1325	c
1326	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1327	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1328	CALL abort_gcm (modname,abort_message,1)
1329	ENDIF
1330	c
1331	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1332	abort_message = 'ISCCP-like outputs may be available for KE
1333	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1334	CALL abort_gcm (modname,abort_message,1)
1335	ENDIF
1336	c
1337	c Initialiser les compteurs:
1338	c
1339	itap = 0
1340	itaprad = 0
1341
1342	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1343	!! Un petit travail à faire ici.
1344	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1345
1346	if (iflag_pbl>1) then
1347	PRINT*, "Using method MELLOR&YAMADA"
1348	endif
1349
1350	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1351	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1352	! dyn3d
1353	! Attention : la version precedente n'etait pas tres propre.
1354	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1355	! pour obtenir le meme resultat.
1356	dtime=pdtphys
1357	radpas = NINT( 86400./dtime/nbapp_rad)
1358	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1359
1360	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1361	cIM begin
1362	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1363	$,ratqs(1,1)
1364	cIM end
1365
1366
1367
1368	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1369	c
1370	C on remet le calendrier a zero
1371	c
1372	IF (raz_date .eq. 1) THEN
1373	itau_phy = 0
1374	ENDIF
1375
1376	cIM cf. AM 081204 BEG
1377	PRINT*,'cycle_diurne3 =',cycle_diurne
1378	cIM cf. AM 081204 END
1379	c
1380	CALL printflag( tabcntr0,radpas,ok_journe,
1381	, ok_instan, ok_region )
1382	c
1383	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1384	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1385	. pdtphys
1386	abort_message='Pas physique n est pas correct '
1387	! call abort_gcm(modname,abort_message,1)
1388	dtime=pdtphys
1389	ENDIF
1390	IF (nlon .NE. klon) THEN
1391	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1392	. klon
1393	abort_message='nlon et klon ne sont pas coherents'
1394	call abort_gcm(modname,abort_message,1)
1395	ENDIF
1396	IF (nlev .NE. klev) THEN
1397	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1398	. klev
1399	abort_message='nlev et klev ne sont pas coherents'
1400	call abort_gcm(modname,abort_message,1)
1401	ENDIF
1402	c
1403	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1404	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1405	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1406	abort_message='Nbre d appels au rayonnement insuffisant'
1407	call abort_gcm(modname,abort_message,1)
1408	ENDIF
1409	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1410	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1411	. ok_cvl
1412	c
1413	cKE43
1414	c Initialisation pour la convection de K.E. (sb):
1415	IF (iflag_con.GE.3) THEN
1416
1417	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1418	WRITE(lunout,*)
1419	. "On va utiliser le melange convectif des traceurs qui"
1420	WRITE(lunout,*)"est calcule dans convect4.3"
1421	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1422
1423	DO i = 1, klon
1424	ema_cbmf(i) = 0.
1425	ema_pcb(i) = 0.
1426	ema_pct(i) = 0.
1427	c ema_workcbmf(i) = 0.
1428	ENDDO
1429	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1430	DO i = 1, klon
1431	ibas_con(i) = 1
1432	itop_con(i) = 1
1433	ENDDO
1434	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1435	c===============================================================================
1436	cCR:04.12.07: initialisations poches froides
1437	c Controle de ALE et ALP pour la fermeture convective (jyg)
1438	if (iflag_wake.eq.1) then
1439	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1440	s ,alp_bl_prescr, ale_bl_prescr)
1441	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1442	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1443	endif
1444
1445	do i = 1,klon
1446	Ale_bl(i)=0.
1447	Alp_bl(i)=0.
1448	enddo
1449
1450	c================================================================================
1451
1452	ENDIF !debut
1453
1454	DO i=1,klon
1455	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1456	ENDDO
1457
1458	c34EK
1459	IF (ok_orodr) THEN
1460
1461	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1462	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1463	! justement quand ok_orodr = false.
1464	! ce rugoro est utilise par la couche limite et fait double emploi
1465	! avec les paramétrisations spécifiques de Francois Lott.
1466	! DO i=1,klon
1467	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1468	! ENDDO
1469	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1470	IF (ok_strato) THEN
1471	CALL SUGWD_strato(klon,klev,paprs,pplay)
1472	ELSE
1473	CALL SUGWD(klon,klev,paprs,pplay)
1474	ENDIF
1475
1476	DO i=1,klon
1477	zuthe(i)=0.
1478	zvthe(i)=0.
1479	if(zstd(i).gt.10.)then
1480	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1481	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1482	endif
1483	ENDDO
1484	ENDIF
1485	c
1486	c
1487	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1488	WRITE(lunout,*)'La frequence de lecture surface est de ',
1489	. lmt_pas
1490	c
1491	cIM 030306 END
1492
1493	capemaxcels = 't_max(X)'
1494	t2mincels = 't_min(X)'
1495	t2maxcels = 't_max(X)'
1496	tinst = 'inst(X)'
1497	tave = 'ave(X)'
1498	cIM cf. AM 081204 BEG
1499	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1500	cIM cf. AM 081204 END
1501	c
1502	c=============================================================
1503	c Initialisation des sorties
1504	c=============================================================
1505
1506	#ifdef CPP_IOIPSL
1507
1508	c$OMP MASTER
1509	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1510	& ctetaSTD,dtime,ok_veget,
1511	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1512	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1513	& read_climoz, new_aod, aerosol_couple
1514	& )
1515	c$OMP END MASTER
1516	c$OMP BARRIER
1517
1518	#ifdef histISCCP
1519	#include "ini_histISCCP.h"
1520	#endif
1521
1522	#ifdef histNMC
1523	#include "ini_histhfNMC.h"
1524	#include "ini_histdayNMC.h"
1525	#include "ini_histmthNMC.h"
1526	#endif
1527
1528	#include "ini_histday_seri.h"
1529
1530	#include "ini_paramLMDZ_phy.h"
1531
1532	#endif
1533
1534	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1535	ecrit_hf2mth = ecrit_mth/ecrit_hf
1536
1537	ecrit_hf = ecrit_hf * un_jour
1538	cIM
1539	IF(ecrit_day.LE.1.) THEN
1540	ecrit_day = ecrit_day * un_jour !en secondes
1541	ENDIF
1542	cIM
1543	ecrit_mth = ecrit_mth * un_jour
1544	ecrit_ins = ecrit_ins * un_jour
1545	ecrit_reg = ecrit_reg * un_jour
1546	ecrit_tra = ecrit_tra * un_jour
1547	ecrit_ISCCP = ecrit_ISCCP * un_jour
1548	ecrit_LES = ecrit_LES * un_jour
1549	c
1550	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1551	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1552	. ecrit_hf2mth
1553	cIM 030306 END
1554
1555
1556	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1557	date0 = jD_ref
1558	WRITE(,) 'physiq date0 : ',date0
1559	c
1560	c
1561	c
1562	c Prescrire l'ozone dans l'atmosphere
1563	c
1564	c
1565	cc DO i = 1, klon
1566	cc DO k = 1, klev
1567	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1568	cc ENDDO
1569	cc ENDDO
1570	c
1571	IF (config_inca /= 'none') THEN
1572	#ifdef INCA
1573	CALL VTe(VTphysiq)
1574	CALL VTb(VTinca)
1575	! iii = MOD(NINT(xjour),360)
1576	! calday = REAL(iii) + jH_cur
1577	calday = REAL(days_elapsed) + jH_cur
1578	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1579
1580	CALL chemini(
1581	$ rg,
1582	$ ra,
1583	$ airephy,
1584	$ rlat,
1585	$ rlon,
1586	$ presnivs,
1587	$ calday,
1588	$ klon,
1589	$ nqtot,
1590	$ pdtphys,
1591	$ annee_ref,
1592	$ day_ref,
1593	$ itau_phy)
1594
1595	CALL VTe(VTinca)
1596	CALL VTb(VTphysiq)
1597	#endif
1598	END IF
1599	c
1600	c
1601	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1602	! Nouvelle initialisation pour le rayonnement RRTM
1603	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1604
1605	call iniradia(klon,klev,paprs(1,1:klev+1))
1606
1607	C$omp single
1608	if (read_climoz >= 1) then
1609	call open_climoz(ncid_climoz, press_climoz)
1610	END IF
1611	C$omp end single
1612	ENDIF
1613	!
1614	! ************** Fin de IF ( debut ) *************
1615	!
1616	!
1617	! Incrementer le compteur de la physique
1618	!
1619	itap = itap + 1
1620	!
1621	! Update fraction of the sub-surfaces (pctsrf) and
1622	! initialize, where a new fraction has appeared, all variables depending
1623	! on the surface fraction.
1624	!
1625	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1626	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1627
1628	! Tendances bidons pour les processus qui n'affectent pas certaines
1629	! variables.
1630	du0(:,:)=0.
1631	dv0(:,:)=0.
1632	dq0(:,:)=0.
1633	dql0(:,:)=0.
1634	c
1635	c Mettre a zero des variables de sortie (pour securite)
1636	c
1637	DO i = 1, klon
1638	d_ps(i) = 0.0
1639	ENDDO
1640	DO k = 1, klev
1641	DO i = 1, klon
1642	d_t(i,k) = 0.0
1643	d_u(i,k) = 0.0
1644	d_v(i,k) = 0.0
1645	ENDDO
1646	ENDDO
1647	DO iq = 1, nqtot
1648	DO k = 1, klev
1649	DO i = 1, klon
1650	d_qx(i,k,iq) = 0.0
1651	ENDDO
1652	ENDDO
1653	ENDDO
1654	da(:,:)=0.
1655	mp(:,:)=0.
1656	phi(:,:,:)=0.
1657	c
1658	c Ne pas affecter les valeurs entrees de u, v, h, et q
1659	c
1660	DO k = 1, klev
1661	DO i = 1, klon
1662	t_seri(i,k) = t(i,k)
1663	u_seri(i,k) = u(i,k)
1664	v_seri(i,k) = v(i,k)
1665	q_seri(i,k) = qx(i,k,ivap)
1666	ql_seri(i,k) = qx(i,k,iliq)
1667	qs_seri(i,k) = 0.
1668	ENDDO
1669	ENDDO
1670	IF (nqtot.GE.3) THEN
1671	DO iq = 3, nqtot
1672	DO k = 1, klev
1673	DO i = 1, klon
1674	tr_seri(i,k,iq-2) = qx(i,k,iq)
1675	ENDDO
1676	ENDDO
1677	ENDDO
1678	ELSE
1679	DO k = 1, klev
1680	DO i = 1, klon
1681	tr_seri(i,k,1) = 0.0
1682	ENDDO
1683	ENDDO
1684	ENDIF
1685	C
1686	DO i = 1, klon
1687	ztsol(i) = 0.
1688	ENDDO
1689	DO nsrf = 1, nbsrf
1690	DO i = 1, klon
1691	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1692	ENDDO
1693	ENDDO
1694	cIM
1695	IF (ip_ebil_phy.ge.1) THEN
1696	ztit='after dynamic'
1697	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1698	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1699	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1700	C Comme les tendances de la physique sont ajoute dans la dynamique,
1701	C on devrait avoir que la variation d'entalpie par la dynamique
1702	C est egale a la variation de la physique au pas de temps precedent.
1703	C Donc la somme de ces 2 variations devrait etre nulle.
1704	call diagphy(airephy,ztit,ip_ebil_phy
1705	e , zero_v, zero_v, zero_v, zero_v, zero_v
1706	e , zero_v, zero_v, zero_v, ztsol
1707	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1708	s , fs_bound, fq_bound )
1709	END IF
1710
1711	c Diagnostiquer la tendance dynamique
1712	c
1713	IF (ancien_ok) THEN
1714	DO k = 1, klev
1715	DO i = 1, klon
1716	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1717	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1718	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1719	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1720	ENDDO
1721	ENDDO
1722	ELSE
1723	DO k = 1, klev
1724	DO i = 1, klon
1725	d_u_dyn(i,k) = 0.0
1726	d_v_dyn(i,k) = 0.0
1727	d_t_dyn(i,k) = 0.0
1728	d_q_dyn(i,k) = 0.0
1729	ENDDO
1730	ENDDO
1731	ancien_ok = .TRUE.
1732	ENDIF
1733	c
1734	c Ajouter le geopotentiel du sol:
1735	c
1736	DO k = 1, klev
1737	DO i = 1, klon
1738	zphi(i,k) = pphi(i,k) + pphis(i)
1739	ENDDO
1740	ENDDO
1741	c
1742	c Verifier les temperatures
1743	c
1744	cIM BEG
1745	IF (check) THEN
1746	amn=MIN(ftsol(1,is_ter),1000.)
1747	amx=MAX(ftsol(1,is_ter),-1000.)
1748	DO i=2, klon
1749	amn=MIN(ftsol(i,is_ter),amn)
1750	amx=MAX(ftsol(i,is_ter),amx)
1751	ENDDO
1752	c
1753	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1754	ENDIF !(check) THEN
1755	cIM END
1756	c
1757	CALL hgardfou(t_seri,ftsol,'debutphy')
1758	c
1759	cIM BEG
1760	IF (check) THEN
1761	amn=MIN(ftsol(1,is_ter),1000.)
1762	amx=MAX(ftsol(1,is_ter),-1000.)
1763	DO i=2, klon
1764	amn=MIN(ftsol(i,is_ter),amn)
1765	amx=MAX(ftsol(i,is_ter),amx)
1766	ENDDO
1767	c
1768	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1769	ENDIF !(check) THEN
1770	cIM END
1771	c
1772	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1773	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1774	c
1775	if (read_climoz >= 1) then
1776	C Ozone from a file
1777	! Update required ozone index:
1778	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1779	$ / ioget_year_len(year_cur) * 360.) + 1
1780	if (ro3i == 361) ro3i = 360
1781	C (This should never occur, except perhaps because of roundup
1782	C error. See documentation.)
1783	if (ro3i /= co3i) then
1784	C Update ozone field:
1785	if (read_climoz == 1) then
1786	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1787	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1788	else
1789	C read_climoz == 2
1790	call regr_pr_av(ncid_climoz,
1791	$ (/"tro3 ", "tro3_daylight"/),
1792	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1793	$ v3=wo)
1794	end if
1795	! Convert from mole fraction of ozone to column density of ozone in a
1796	! cell, in kDU:
1797	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1798	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1799	C (By regridding ozone values for LMDZ only once every 360th of
1800	C year, we have already neglected the variation of pressure in one
1801	C 360th of year. So do not recompute "wo" at each time step even if
1802	C "zmasse" changes a little.)
1803	co3i = ro3i
1804	end if
1805	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1806	C Once per day, update ozone from Royer:
1807	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1808	ENDIF
1809	c
1810	c Re-evaporer l'eau liquide nuageuse
1811	c
1812	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1813	DO i = 1, klon
1814	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1815	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1816	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1817	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1818	zb = MAX(0.0,ql_seri(i,k))
1819	za = - MAX(0.0,ql_seri(i,k))
1820	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1821	t_seri(i,k) = t_seri(i,k) + za
1822	q_seri(i,k) = q_seri(i,k) + zb
1823	ql_seri(i,k) = 0.0
1824	d_t_eva(i,k) = za
1825	d_q_eva(i,k) = zb
1826	ENDDO
1827	ENDDO
1828	cIM
1829	IF (ip_ebil_phy.ge.2) THEN
1830	ztit='after reevap'
1831	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1832	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1833	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1834	call diagphy(airephy,ztit,ip_ebil_phy
1835	e , zero_v, zero_v, zero_v, zero_v, zero_v
1836	e , zero_v, zero_v, zero_v, ztsol
1837	e , d_h_vcol, d_qt, d_ec
1838	s , fs_bound, fq_bound )
1839	C
1840	END IF
1841
1842	c
1843	c=========================================================================
1844	! Calculs de l'orbite.
1845	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1846	! doit donc etre placé avant radlwsw et pbl_surface
1847
1848	! calcul selon la routine utilisee pour les planetes
1849	if (new_orbit) then
1850	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1851	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1852	! day_since_equinox = (jD_cur) - jD_eq
1853	call solarlong(day_since_equinox, zlongi, dist)
1854	else
1855	! calcul selon la routine utilisee pour l'AR4
1856	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1857	! solarlong0
1858	if (solarlong0<-999.) then
1859	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1860	else
1861	zlongi=solarlong0 ! longitude solaire vraie
1862	dist=1. ! distance au soleil / moyenne
1863	endif
1864	endif
1865	if(prt_level.ge.1) &
1866	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1867
1868	! Avec ou sans cycle diurne
1869	IF (cycle_diurne) THEN
1870	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1871	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1872	ELSE
1873	CALL angle(zlongi, rlat, fract, rmu0)
1874	ENDIF
1875
1876	if (mydebug) then
1877	call writefield_phy('u_seri',u_seri,llm)
1878	call writefield_phy('v_seri',v_seri,llm)
1879	call writefield_phy('t_seri',t_seri,llm)
1880	call writefield_phy('q_seri',q_seri,llm)
1881	endif
1882
1883	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1884	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1885	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1886	c turbulent du couche limit.
1887	c
1888	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1889	c ecriture des fihiers hist_XXXX.nc, ces sont :
1890	c qsol, zq2m, s_pblh, s_lcl,
1891	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1892	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1893	c zxrugs, zu10m, zv10m, fder,
1894	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1895	c frugs, agesno, fsollw, fsolsw,
1896	c d_ts, fevap, fluxlat, t2m,
1897	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1898	c
1899	c Certains ne sont pas utiliser du tout :
1900	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1901	c
1902
1903	CALL pbl_surface(
1904	e dtime, date0, itap, days_elapsed+1,
1905	e debut, lafin,
1906	e rlon, rlat, rugoro, rmu0,
1907	e rain_fall, snow_fall, solsw, sollw,
1908	e t_seri, q_seri, u_seri, v_seri,
1909	e pplay, paprs, pctsrf,
1910	+ ftsol, falb1, falb2, u10m, v10m,
1911	s sollwdown, cdragh, cdragm, u1, v1,
1912	s albsol1, albsol2, sens, evap,
1913	s zxtsol, zxfluxlat, zt2m, qsat2m,
1914	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1915	s coefh, slab_wfbils,
1916	d qsol, zq2m, s_pblh, s_lcl,
1917	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1918	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1919	d zxrugs, zu10m, zv10m, fder,
1920	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1921	d frugs, agesno, fsollw, fsolsw,
1922	d d_ts, fevap, fluxlat, t2m,
1923	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1924	- dsens, devap, zxsnow,
1925	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1926
1927
1928	!-----------------------------------------------------------------------------------------
1929	! ajout des tendances de la diffusion turbulente
1930	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1931	!-----------------------------------------------------------------------------------------
1932
1933	if (mydebug) then
1934	call writefield_phy('u_seri',u_seri,llm)
1935	call writefield_phy('v_seri',v_seri,llm)
1936	call writefield_phy('t_seri',t_seri,llm)
1937	call writefield_phy('q_seri',q_seri,llm)
1938	endif
1939
1940
1941	IF (ip_ebil_phy.ge.2) THEN
1942	ztit='after surface_main'
1943	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1944	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1945	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1946	call diagphy(airephy,ztit,ip_ebil_phy
1947	e , zero_v, zero_v, zero_v, zero_v, sens
1948	e , evap , zero_v, zero_v, ztsol
1949	e , d_h_vcol, d_qt, d_ec
1950	s , fs_bound, fq_bound )
1951	END IF
1952
1953	c =================================================================== c
1954	c Calcul de Qsat
1955
1956	DO k = 1, klev
1957	DO i = 1, klon
1958	zx_t = t_seri(i,k)
1959	IF (thermcep) THEN
1960	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1961	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1962	zx_qs = MIN(0.5,zx_qs)
1963	zcor = 1./(1.-retv*zx_qs)
1964	zx_qs = zx_qs*zcor
1965	ELSE
1966	IF (zx_t.LT.t_coup) THEN
1967	zx_qs = qsats(zx_t)/pplay(i,k)
1968	ELSE
1969	zx_qs = qsatl(zx_t)/pplay(i,k)
1970	ENDIF
1971	ENDIF
1972	zqsat(i,k)=zx_qs
1973	ENDDO
1974	ENDDO
1975
1976	if (prt_level.ge.1) then
1977	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1978	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1979	endif
1980	c
1981	c Appeler la convection (au choix)
1982	c
1983	DO k = 1, klev
1984	DO i = 1, klon
1985	conv_q(i,k) = d_q_dyn(i,k)
1986	. + d_q_vdf(i,k)/dtime
1987	conv_t(i,k) = d_t_dyn(i,k)
1988	. + d_t_vdf(i,k)/dtime
1989	ENDDO
1990	ENDDO
1991	IF (check) THEN
1992	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1993	WRITE(lunout,*) "avantcon=", za
1994	ENDIF
1995	zx_ajustq = .FALSE.
1996	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1997	IF (zx_ajustq) THEN
1998	DO i = 1, klon
1999	z_avant(i) = 0.0
2000	ENDDO
2001	DO k = 1, klev
2002	DO i = 1, klon
2003	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2004	. *(paprs(i,k)-paprs(i,k+1))/RG
2005	ENDDO
2006	ENDDO
2007	ENDIF
2008
2009	c Calcule de vitesse verticale a partir de flux de masse verticale
2010	DO k = 1, klev
2011	DO i = 1, klon
2012	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2013	END DO
2014	END DO
2015	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
2016	$ omega(igout, :)
2017
2018	IF (iflag_con.EQ.1) THEN
2019	abort_message ='reactiver le call conlmd dans physiq.F'
2020	CALL abort_gcm (modname,abort_message,1)
2021	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2022	c . d_t_con, d_q_con,
2023	c . rain_con, snow_con, ibas_con, itop_con)
2024	ELSE IF (iflag_con.EQ.2) THEN
2025	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2026	e conv_t, conv_q, -evap, omega,
2027	s d_t_con, d_q_con, rain_con, snow_con,
2028	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2029	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2030	d_u_con = 0.
2031	d_v_con = 0.
2032
2033	WHERE (rain_con < 0.) rain_con = 0.
2034	WHERE (snow_con < 0.) snow_con = 0.
2035	DO i = 1, klon
2036	ibas_con(i) = klev+1 - kcbot(i)
2037	itop_con(i) = klev+1 - kctop(i)
2038	ENDDO
2039	ELSE IF (iflag_con.GE.3) THEN
2040	c nb of tracers for the KE convection:
2041	c MAF la partie traceurs est faite dans phytrac
2042	c on met ntra=1 pour limiter les appels mais on peut
2043	c supprimer les calculs / ftra.
2044	ntra = 1
2045
2046	c=====================================================================================
2047	cajout pour la parametrisation des poches froides:
2048	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2049	do k=1,klev
2050	do i=1,klon
2051	if (iflag_wake.eq.1) then
2052	t_wake(i,k) = t_seri(i,k)
2053	. +(1-wake_s(i))*wake_deltat(i,k)
2054	q_wake(i,k) = q_seri(i,k)
2055	. +(1-wake_s(i))*wake_deltaq(i,k)
2056	t_undi(i,k) = t_seri(i,k)
2057	. -wake_s(i)*wake_deltat(i,k)
2058	q_undi(i,k) = q_seri(i,k)
2059	. -wake_s(i)*wake_deltaq(i,k)
2060	else
2061	t_wake(i,k) = t_seri(i,k)
2062	q_wake(i,k) = q_seri(i,k)
2063	t_undi(i,k) = t_seri(i,k)
2064	q_undi(i,k) = q_seri(i,k)
2065	endif
2066	enddo
2067	enddo
2068
2069	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2070	cc-- pour le soulevement des particules dans le modele convectif
2071	c
2072	do i = 1,klon
2073	ALE(i) = 0.
2074	ALP(i) = 0.
2075	enddo
2076	c
2077	ccalcul de ale_wake et alp_wake
2078	if (iflag_wake.eq.1) then
2079	if (itap .le. it_wape_prescr) then
2080	do i = 1,klon
2081	ale_wake(i) = wape_prescr
2082	alp_wake(i) = fip_prescr
2083	enddo
2084	else
2085	do i = 1,klon
2086	cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2087	ccc ale_wake(i) = 0.5wake_cstar(i)*2
2088	ale_wake(i) = wake_pe(i)
2089	alp_wake(i) = wake_fip(i)
2090	enddo
2091	endif
2092	else
2093	do i = 1,klon
2094	ale_wake(i) = 0.
2095	alp_wake(i) = 0.
2096	enddo
2097	endif
2098	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2099	cdans le thermique sinon
2100	if (iflag_coupl.eq.0) then
2101	if (debut.and.prt_level.gt.9)
2102	$ WRITE(lunout,*)'ALE et ALP imposes'
2103	do i = 1,klon
2104	con ne couple que ale
2105	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2106	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2107	con ne couple que alp
2108	c ALP(i) = alp_wake(i) + Alp_bl(i)
2109	ALP(i) = alp_wake(i) + alp_bl_prescr
2110	enddo
2111	else
2112	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2113	! do i = 1,klon
2114	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2115	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2116	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2117	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2118	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2119	! enddo
2120
2121	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2122	! Modif FH 2010/04/27. Sans doute temporaire.
2123	! Deux options pour le alp_offset : constant si >Ã 0 ou proportionnel Ãa
2124	! w si <0
2125	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2126	do i = 1,klon
2127	ALE(i) = max(ale_wake(i),Ale_bl(i))
2128	if (alp_offset>=0.) then
2129	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2130	else
2131	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2132	if (alp(i)<0.) then
2133	print*,'ALP ',alp(i),alp_wake(i)
2134	s ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2135	endif
2136	endif
2137	enddo
2138	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2139
2140	endif
2141	do i=1,klon
2142	if (alp(i)>alp_max) then
2143	IF(prt_level>9)WRITE(lunout,*) &
2144	& 'WARNING SUPER ALP (seuil=',alp_max,
2145	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2146	alp(i)=alp_max
2147	endif
2148	if (ale(i)>ale_max) then
2149	IF(prt_level>9)WRITE(lunout,*) &
2150	& 'WARNING SUPER ALE (seuil=',ale_max,
2151	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2152	ale(i)=ale_max
2153	endif
2154	enddo
2155
2156	cfin calcul ale et alp
2157	c=================================================================================================
2158
2159
2160	c sb, oct02:
2161	c Schema de convection modularise et vectorise:
2162	c (driver commun aux versions 3 et 4)
2163	c
2164	IF (ok_cvl) THEN ! new driver for convectL
2165
2166	CALL concvl (iflag_con,iflag_clos,
2167	. dtime,paprs,pplay,t_undi,q_undi,
2168	. t_wake,q_wake,wake_s,
2169	. u_seri,v_seri,tr_seri,nbtr,
2170	. ALE,ALP,
2171	. ema_work1,ema_work2,
2172	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2173	. rain_con, snow_con, ibas_con, itop_con, sigd,
2174	. ema_cbmf,upwd,dnwd,dnwd0,
2175	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2176	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2177	. pmflxr,pmflxs,da,phi,mp,
2178	. ftd,fqd,lalim_conv,wght_th)
2179
2180	cIM begin
2181	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2182	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2183	cIM end
2184	cIM cf. FH
2185	clwcon0=qcondc
2186	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2187
2188	do i = 1, klon
2189	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2190	enddo
2191
2192	ELSE ! ok_cvl
2193
2194	c MAF conema3 ne contient pas les traceurs
2195	CALL conema3 (dtime,
2196	. paprs,pplay,t_seri,q_seri,
2197	. u_seri,v_seri,tr_seri,ntra,
2198	. ema_work1,ema_work2,
2199	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2200	. rain_con, snow_con, ibas_con, itop_con,
2201	. upwd,dnwd,dnwd0,bas,top,
2202	. Ma,cape,tvp,rflag,
2203	. pbase
2204	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2205	. ,clwcon0)
2206
2207	ENDIF ! ok_cvl
2208
2209	c
2210	c Correction precip
2211	rain_con = rain_con * cvl_corr
2212	snow_con = snow_con * cvl_corr
2213	c
2214
2215	IF (.NOT. ok_gust) THEN
2216	do i = 1, klon
2217	wd(i)=0.0
2218	enddo
2219	ENDIF
2220
2221	c =================================================================== c
2222	c Calcul des proprietes des nuages convectifs
2223	c
2224
2225	c calcul des proprietes des nuages convectifs
2226	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2227	call clouds_gno
2228	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2229
2230	c =================================================================== c
2231
2232	DO i = 1, klon
2233	itop_con(i) = min(max(itop_con(i),1),klev)
2234	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2235	ENDDO
2236
2237	DO i = 1, klon
2238	ema_pcb(i) = paprs(i,ibas_con(i))
2239	ENDDO
2240	DO i = 1, klon
2241	! L'idicage de itop_con peut cacher un pb potentiel
2242	! FH sous la dictee de JYG, CR
2243	ema_pct(i) = paprs(i,itop_con(i)+1)
2244
2245	if (itop_con(i).gt.klev-3) then
2246	if(prt_level >= 9) then
2247	write(lunout,*)'La convection monte trop haut '
2248	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2249	endif
2250	endif
2251	ENDDO
2252	ELSE IF (iflag_con.eq.0) THEN
2253	write(lunout,*) 'On n appelle pas la convection'
2254	clwcon0=0.
2255	rnebcon0=0.
2256	d_t_con=0.
2257	d_q_con=0.
2258	d_u_con=0.
2259	d_v_con=0.
2260	rain_con=0.
2261	snow_con=0.
2262	bas=1
2263	top=1
2264	ELSE
2265	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2266	CALL abort
2267	ENDIF
2268
2269	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2270	c . d_u_con, d_v_con)
2271
2272	!-----------------------------------------------------------------------------------------
2273	! ajout des tendances de la diffusion turbulente
2274	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2275	!-----------------------------------------------------------------------------------------
2276
2277	if (mydebug) then
2278	call writefield_phy('u_seri',u_seri,llm)
2279	call writefield_phy('v_seri',v_seri,llm)
2280	call writefield_phy('t_seri',t_seri,llm)
2281	call writefield_phy('q_seri',q_seri,llm)
2282	endif
2283
2284	cIM
2285	IF (ip_ebil_phy.ge.2) THEN
2286	ztit='after convect'
2287	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2288	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2289	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2290	call diagphy(airephy,ztit,ip_ebil_phy
2291	e , zero_v, zero_v, zero_v, zero_v, zero_v
2292	e , zero_v, rain_con, snow_con, ztsol
2293	e , d_h_vcol, d_qt, d_ec
2294	s , fs_bound, fq_bound )
2295	END IF
2296	C
2297	IF (check) THEN
2298	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2299	WRITE(lunout,*)"aprescon=", za
2300	zx_t = 0.0
2301	za = 0.0
2302	DO i = 1, klon
2303	za = za + airephy(i)/REAL(klon)
2304	zx_t = zx_t + (rain_con(i)+
2305	. snow_con(i))*airephy(i)/REAL(klon)
2306	ENDDO
2307	zx_t = zx_t/za*dtime
2308	WRITE(lunout,*)"Precip=", zx_t
2309	ENDIF
2310	IF (zx_ajustq) THEN
2311	DO i = 1, klon
2312	z_apres(i) = 0.0
2313	ENDDO
2314	DO k = 1, klev
2315	DO i = 1, klon
2316	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2317	. *(paprs(i,k)-paprs(i,k+1))/RG
2318	ENDDO
2319	ENDDO
2320	DO i = 1, klon
2321	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2322	. /z_apres(i)
2323	ENDDO
2324	DO k = 1, klev
2325	DO i = 1, klon
2326	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2327	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2328	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2329	ENDIF
2330	ENDDO
2331	ENDDO
2332	ENDIF
2333	zx_ajustq=.FALSE.
2334
2335	c
2336	c=============================================================================
2337	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2338	cpour la couche limite diffuse pour l instant
2339	c
2340	if (iflag_wake.eq.1) then
2341	DO k=1,klev
2342	DO i=1,klon
2343	dt_dwn(i,k) = ftd(i,k)
2344	wdt_PBL(i,k) = 0.
2345	dq_dwn(i,k) = fqd(i,k)
2346	wdq_PBL(i,k) = 0.
2347	M_dwn(i,k) = dnwd0(i,k)
2348	M_up(i,k) = upwd(i,k)
2349	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2350	udt_PBL(i,k) = 0.
2351	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2352	udq_PBL(i,k) = 0.
2353	ENDDO
2354	ENDDO
2355	c
2356	ccalcul caracteristiques de la poche froide
2357	call calWAKE (paprs,pplay,dtime
2358	: ,t_seri,q_seri,omega
2359	: ,dt_dwn,dq_dwn,M_dwn,M_up
2360	: ,dt_a,dq_a,sigd
2361	: ,wdt_PBL,wdq_PBL
2362	: ,udt_PBL,udq_PBL
2363	o ,wake_deltat,wake_deltaq,wake_dth
2364	o ,wake_h,wake_s,wake_dens
2365	o ,wake_pe,wake_fip,wake_gfl
2366	o ,dt_wake,dq_wake
2367	o ,wake_k, t_undi,q_undi
2368	o ,wake_omgbdth,wake_dp_omgb
2369	o ,wake_dtKE,wake_dqKE
2370	o ,wake_dtPBL,wake_dqPBL
2371	o ,wake_omg,wake_dp_deltomg
2372	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2373	o ,wake_ddeltat,wake_ddeltaq)
2374	c
2375	!-----------------------------------------------------------------------------------------
2376	! ajout des tendances des poches froides
2377	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2378	! coherent avec les autres d_t_...
2379	d_t_wake(:,:)=dt_wake(:,:)*dtime
2380	d_q_wake(:,:)=dq_wake(:,:)*dtime
2381	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2382	!-----------------------------------------------------------------------------------------
2383
2384	endif
2385	c
2386	c===================================================================
2387	cJYG
2388	IF (ip_ebil_phy.ge.2) THEN
2389	ztit='after wake'
2390	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2391	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2392	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2393	call diagphy(airephy,ztit,ip_ebil_phy
2394	e , zero_v, zero_v, zero_v, zero_v, zero_v
2395	e , zero_v, zero_v, zero_v, ztsol
2396	e , d_h_vcol, d_qt, d_ec
2397	s , fs_bound, fq_bound )
2398	END IF
2399
2400	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2401	c
2402	c===================================================================
2403	c Convection seche (thermiques ou ajustement)
2404	c===================================================================
2405	c
2406	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2407	s ,seuil_inversion,weak_inversion,dthmin)
2408
2409
2410
2411	d_t_ajsb(:,:)=0.
2412	d_q_ajsb(:,:)=0.
2413	d_t_ajs(:,:)=0.
2414	d_u_ajs(:,:)=0.
2415	d_v_ajs(:,:)=0.
2416	d_q_ajs(:,:)=0.
2417	clwcon0th(:,:)=0.
2418	c
2419	c fm_therm(:,:)=0.
2420	c entr_therm(:,:)=0.
2421	c detr_therm(:,:)=0.
2422	c
2423	IF(prt_level>9)WRITE(lunout,*)
2424	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2425	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2426	if(iflag_thermals.lt.0) then
2427	c Rien
2428	c ====
2429	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2430
2431
2432	else
2433
2434	c Thermiques
2435	c ==========
2436	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2437	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2438
2439
2440	if (iflag_thermals.gt.1) then
2441	call calltherm(pdtphys
2442	s ,pplay,paprs,pphi,weak_inversion
2443	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2444	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2445	s ,fm_therm,entr_therm,detr_therm
2446	s ,zqasc,clwcon0th,lmax_th,ratqscth
2447	s ,ratqsdiff,zqsatth
2448	con rajoute ale et alp, et les caracteristiques de la couche alim
2449	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
2450	s ,ztv,zpspsk,ztla,zthl)
2451
2452	! ----------------------------------------------------------------------
2453	! Transport de la TKE par les panaches thermiques.
2454	! FH : 2010/02/01
2455	if (iflag_pbl.eq.10) then
2456	call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2457	s rg,paprs,pbl_tke)
2458	endif
2459	! ----------------------------------------------------------------------
2460
2461	endif
2462
2463
2464
2465	c Ajustement sec
2466	c ==============
2467
2468	! Dans le cas oÃ¹ on active les thermiques, on fait partir l'ajustement
2469	! a partir du sommet des thermiques.
2470	! Dans le cas contraire, on demarre au niveau 1.
2471
2472	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2473
2474	if(iflag_thermals.eq.0) then
2475	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2476	limbas(:)=1
2477	else
2478	limbas(:)=lmax_th(:)
2479	endif
2480
2481	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2482	! pour des test de convergence numerique.
2483	! Le nouveau ajsec est a priori mieux, meme pour le cas
2484	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2485	! non nulles numeriquement pour des mailles non concernees.
2486
2487	if (iflag_thermals.eq.0) then
2488	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2489	s , d_t_ajsb, d_q_ajsb)
2490	else
2491	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2492	s , d_t_ajsb, d_q_ajsb)
2493	endif
2494
2495	!-----------------------------------------------------------------------------------------
2496	! ajout des tendances de l'ajustement sec ou des thermiques
2497	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2498	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2499	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2500
2501	!-----------------------------------------------------------------------------------------
2502
2503	endif
2504
2505	endif
2506	c
2507	c===================================================================
2508	cIM
2509	IF (ip_ebil_phy.ge.2) THEN
2510	ztit='after dry_adjust'
2511	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2512	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2513	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2514	call diagphy(airephy,ztit,ip_ebil_phy
2515	e , zero_v, zero_v, zero_v, zero_v, zero_v
2516	e , zero_v, zero_v, zero_v, ztsol
2517	e , d_h_vcol, d_qt, d_ec
2518	s , fs_bound, fq_bound )
2519	END IF
2520
2521
2522	c-------------------------------------------------------------------------
2523	c Caclul des ratqs
2524	c-------------------------------------------------------------------------
2525
2526	c print*,'calcul des ratqs'
2527	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2528	c ----------------
2529	c on ecrase le tableau ratqsc calcule par clouds_gno
2530	if (iflag_cldcon.eq.1) then
2531	do k=1,klev
2532	do i=1,klon
2533	if(ptconv(i,k)) then
2534	ratqsc(i,k)=ratqsbas
2535	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2536	else
2537	ratqsc(i,k)=0.
2538	endif
2539	enddo
2540	enddo
2541
2542	c-----------------------------------------------------------------------
2543	c par nversion de la fonction log normale
2544	c-----------------------------------------------------------------------
2545	else if (iflag_cldcon.eq.4) then
2546	ptconvth(:,:)=.false.
2547	ratqsc(:,:)=0.
2548	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2549	call clouds_gno
2550	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2551	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2552
2553	endif
2554
2555	c ratqs stables
2556	c -------------
2557
2558	if (iflag_ratqs.eq.0) then
2559
2560	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2561	do k=1,klev
2562	do i=1, klon
2563	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2564	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2565	enddo
2566	enddo
2567
2568	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2569	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2570	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2571	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2572	! Il s'agit de differents tests dans la phase de reglage du modele
2573	! avec thermiques.
2574
2575	else if (iflag_ratqs.eq.1) then
2576
2577	do k=1,klev
2578	do i=1, klon
2579	if (pplay(i,k).ge.60000.) then
2580	ratqss(i,k)=ratqsbas
2581	else if ((pplay(i,k).ge.30000.).and.
2582	s (pplay(i,k).lt.60000.)) then
2583	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2584	s (60000.-pplay(i,k))/(60000.-30000.)
2585	else
2586	ratqss(i,k)=ratqshaut
2587	endif
2588	enddo
2589	enddo
2590
2591	else
2592
2593	do k=1,klev
2594	do i=1, klon
2595	if (pplay(i,k).ge.60000.) then
2596	ratqss(i,k)=ratqsbas
2597	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2598	else if ((pplay(i,k).ge.30000.).and.
2599	s (pplay(i,k).lt.60000.)) then
2600	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2601	s (60000.-pplay(i,k))/(60000.-30000.)
2602	else
2603	ratqss(i,k)=ratqshaut
2604	endif
2605	enddo
2606	enddo
2607	endif
2608
2609
2610
2611
2612	c ratqs final
2613	c -----------
2614
2615	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2616	s .or.iflag_cldcon.eq.4) then
2617
2618	! On ajoute une constante au ratqsc*2 pour tenir compte de
2619	! fluctuations turbulentes de petite echelle
2620
2621	do k=1,klev
2622	do i=1,klon
2623	if ((fm_therm(i,k).gt.1.e-10)) then
2624	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2625	endif
2626	enddo
2627	enddo
2628
2629	! les ratqs sont une combinaison de ratqss et ratqsc
2630	if(prt_level.ge.9)
2631	$ write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2632
2633	if (tau_ratqs>1.e-10) then
2634	facteur=exp(-pdtphys/tau_ratqs)
2635	else
2636	facteur=0.
2637	endif
2638	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2639	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2640	! FH 22/09/2009
2641	! La ligne ci-dessous faisait osciller le modele et donnait une solution
2642	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
2643	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2644	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2645	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2646	else
2647	! on ne prend que le ratqs stable pour fisrtilp
2648	ratqs(:,:)=ratqss(:,:)
2649	endif
2650
2651
2652	c
2653	c Appeler le processus de condensation a grande echelle
2654	c et le processus de precipitation
2655	c-------------------------------------------------------------------------
2656	CALL fisrtilp(dtime,paprs,pplay,
2657	. t_seri, q_seri,ptconv,ratqs,
2658	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2659	. rain_lsc, snow_lsc,
2660	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2661	. frac_impa, frac_nucl,
2662	. prfl, psfl, rhcl,
2663	. zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
2664
2665	WHERE (rain_lsc < 0) rain_lsc = 0.
2666	WHERE (snow_lsc < 0) snow_lsc = 0.
2667	!-----------------------------------------------------------------------------------------
2668	! ajout des tendances de la diffusion turbulente
2669	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2670	!-----------------------------------------------------------------------------------------
2671	DO k = 1, klev
2672	DO i = 1, klon
2673	cldfra(i,k) = rneb(i,k)
2674	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2675	ENDDO
2676	ENDDO
2677	IF (check) THEN
2678	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2679	WRITE(lunout,*)"apresilp=", za
2680	zx_t = 0.0
2681	za = 0.0
2682	DO i = 1, klon
2683	za = za + airephy(i)/REAL(klon)
2684	zx_t = zx_t + (rain_lsc(i)
2685	. + snow_lsc(i))*airephy(i)/REAL(klon)
2686	ENDDO
2687	zx_t = zx_t/za*dtime
2688	WRITE(lunout,*)"Precip=", zx_t
2689	ENDIF
2690	cIM
2691	IF (ip_ebil_phy.ge.2) THEN
2692	ztit='after fisrt'
2693	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2694	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2695	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2696	call diagphy(airephy,ztit,ip_ebil_phy
2697	e , zero_v, zero_v, zero_v, zero_v, zero_v
2698	e , zero_v, rain_lsc, snow_lsc, ztsol
2699	e , d_h_vcol, d_qt, d_ec
2700	s , fs_bound, fq_bound )
2701	END IF
2702
2703	if (mydebug) then
2704	call writefield_phy('u_seri',u_seri,llm)
2705	call writefield_phy('v_seri',v_seri,llm)
2706	call writefield_phy('t_seri',t_seri,llm)
2707	call writefield_phy('q_seri',q_seri,llm)
2708	endif
2709
2710	c
2711	c-------------------------------------------------------------------
2712	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2713	c-------------------------------------------------------------------
2714
2715	c 1. NUAGES CONVECTIFS
2716	c
2717	cIM cf FH
2718	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2719	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2720	snow_tiedtke=0.
2721	c print*,'avant calcul de la pseudo precip '
2722	c print*,'iflag_cldcon',iflag_cldcon
2723	if (iflag_cldcon.eq.-1) then
2724	rain_tiedtke=rain_con
2725	else
2726	c print*,'calcul de la pseudo precip '
2727	rain_tiedtke=0.
2728	c print*,'calcul de la pseudo precip 0'
2729	do k=1,klev
2730	do i=1,klon
2731	if (d_q_con(i,k).lt.0.) then
2732	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2733	s *(paprs(i,k)-paprs(i,k+1))/rg
2734	endif
2735	enddo
2736	enddo
2737	endif
2738	c
2739	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2740	c
2741
2742	c Nuages diagnostiques pour Tiedtke
2743	CALL diagcld1(paprs,pplay,
2744	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2745	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2746	. diafra,dialiq)
2747	DO k = 1, klev
2748	DO i = 1, klon
2749	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2750	cldliq(i,k) = dialiq(i,k)
2751	cldfra(i,k) = diafra(i,k)
2752	ENDIF
2753	ENDDO
2754	ENDDO
2755
2756	ELSE IF (iflag_cldcon.ge.3) THEN
2757	c On prend pour les nuages convectifs le max du calcul de la
2758	c convection et du calcul du pas de temps precedent diminue d'un facteur
2759	c facttemps
2760	facteur = pdtphys *facttemps
2761	do k=1,klev
2762	do i=1,klon
2763	rnebcon(i,k)=rnebcon(i,k)*facteur
2764	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2765	s then
2766	rnebcon(i,k)=rnebcon0(i,k)
2767	clwcon(i,k)=clwcon0(i,k)
2768	endif
2769	enddo
2770	enddo
2771
2772	c
2773	cjq - introduce the aerosol direct and first indirect radiative forcings
2774	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2775	IF (ok_ade.OR.ok_aie) THEN
2776	IF (.NOT. aerosol_couple)
2777	& CALL readaerosol_optic(
2778	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2779	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2780	& mass_solu_aero, mass_solu_aero_pi,
2781	& tau_aero, piz_aero, cg_aero,
2782	& tausum_aero, tau3d_aero)
2783	ELSE
2784	cIM 170310 BEG
2785	tausum_aero(:,:,:) = 0.
2786	cIM 170310 END
2787	tau_aero(:,:,:,:) = 0.
2788	piz_aero(:,:,:,:) = 0.
2789	cg_aero(:,:,:,:) = 0.
2790	ENDIF
2791
2792	cIM calcul nuages par le simulateur ISCCP
2793	c
2794	#ifdef histISCCP
2795	IF (ok_isccp) THEN
2796	c
2797	cIM lecture invtau, tautab des fichiers formattes
2798	c
2799	IF (debut) THEN
2800	c$OMP MASTER
2801	c
2802	open(99,file='tautab.formatted', FORM='FORMATTED')
2803	read(99,'(f30.20)') tautab_omp
2804	close(99)
2805	c
2806	open(99,file='invtau.formatted',form='FORMATTED')
2807	read(99,'(i10)') invtau_omp
2808
2809	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2810	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2811
2812	close(99)
2813	c$OMP END MASTER
2814	c$OMP BARRIER
2815	tautab=tautab_omp
2816	invtau=invtau_omp
2817	c
2818	ENDIF !debut
2819	c
2820	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2821	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2822	#include "calcul_simulISCCP.h"
2823	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2824	ENDIF !ok_isccp
2825	#endif
2826
2827	c On prend la somme des fractions nuageuses et des contenus en eau
2828
2829	if (iflag_cldcon>=5) then
2830
2831	! Si on est sur un point touche par la convection profonde et pas
2832	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
2833	! de la convection profonde.
2834
2835	print*,'TEST SCHEMA DE NUAGES '
2836	do k=1,klev
2837	do i=1,klon
2838	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
2839	cldfra(i,k)=rnebcon(i,k)
2840	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
2841	endif
2842	enddo
2843	enddo
2844	else
2845	! Ancienne version
2846	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2847	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2848	endif
2849
2850	ENDIF
2851	c
2852	c 2. NUAGES STARTIFORMES
2853	c
2854	IF (ok_stratus) THEN
2855	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2856	DO k = 1, klev
2857	DO i = 1, klon
2858	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2859	cldliq(i,k) = dialiq(i,k)
2860	cldfra(i,k) = diafra(i,k)
2861	ENDIF
2862	ENDDO
2863	ENDDO
2864	ENDIF
2865	c
2866	c Precipitation totale
2867	c
2868	DO i = 1, klon
2869	rain_fall(i) = rain_con(i) + rain_lsc(i)
2870	snow_fall(i) = snow_con(i) + snow_lsc(i)
2871	ENDDO
2872	cIM
2873	IF (ip_ebil_phy.ge.2) THEN
2874	ztit="after diagcld"
2875	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2876	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2877	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2878	call diagphy(airephy,ztit,ip_ebil_phy
2879	e , zero_v, zero_v, zero_v, zero_v, zero_v
2880	e , zero_v, zero_v, zero_v, ztsol
2881	e , d_h_vcol, d_qt, d_ec
2882	s , fs_bound, fq_bound )
2883	END IF
2884	c
2885	c Calculer l'humidite relative pour diagnostique
2886	c
2887	DO k = 1, klev
2888	DO i = 1, klon
2889	zx_t = t_seri(i,k)
2890	IF (thermcep) THEN
2891	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2892	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2893	zx_qs = MIN(0.5,zx_qs)
2894	zcor = 1./(1.-retv*zx_qs)
2895	zx_qs = zx_qs*zcor
2896	ELSE
2897	IF (zx_t.LT.t_coup) THEN
2898	zx_qs = qsats(zx_t)/pplay(i,k)
2899	ELSE
2900	zx_qs = qsatl(zx_t)/pplay(i,k)
2901	ENDIF
2902	ENDIF
2903	zx_rh(i,k) = q_seri(i,k)/zx_qs
2904	zqsat(i,k)=zx_qs
2905	ENDDO
2906	ENDDO
2907
2908	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2909	c equivalente a 2m (tpote) pour diagnostique
2910	c
2911	DO i = 1, klon
2912	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2913	IF (thermcep) THEN
2914	IF(zt2m(i).LT.RTT) then
2915	Lheat=RLSTT
2916	ELSE
2917	Lheat=RLVTT
2918	ENDIF
2919	ELSE
2920	IF (zt2m(i).LT.RTT) THEN
2921	Lheat=RLSTT
2922	ELSE
2923	Lheat=RLVTT
2924	ENDIF
2925	ENDIF
2926	tpote(i) = tpot(i)*
2927	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2928	ENDDO
2929
2930	IF (config_inca /= 'none') THEN
2931	#ifdef INCA
2932	CALL VTe(VTphysiq)
2933	CALL VTb(VTinca)
2934	calday = REAL(days_elapsed + 1) + jH_cur
2935
2936	call chemtime(itap+itau_phy-1, date0, dtime)
2937	IF (config_inca == 'aero') THEN
2938	CALL AEROSOL_METEO_CALC(
2939	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2940	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2941	END IF
2942
2943	zxsnow_dummy(:) = 0.0
2944
2945	CALL chemhook_begin (calday,
2946	$ days_elapsed+1,
2947	$ jH_cur,
2948	$ pctsrf(1,1),
2949	$ rlat,
2950	$ rlon,
2951	$ airephy,
2952	$ paprs,
2953	$ pplay,
2954	$ coefh,
2955	$ pphi,
2956	$ t_seri,
2957	$ u,
2958	$ v,
2959	$ wo(:, :, 1),
2960	$ q_seri,
2961	$ zxtsol,
2962	$ zxsnow_dummy,
2963	$ solsw,
2964	$ albsol1,
2965	$ rain_fall,
2966	$ snow_fall,
2967	$ itop_con,
2968	$ ibas_con,
2969	$ cldfra,
2970	$ iim,
2971	$ jjm,
2972	$ tr_seri,
2973	$ ftsol,
2974	$ paprs,
2975	$ cdragh,
2976	$ cdragm,
2977	$ pctsrf,
2978	$ pdtphys,
2979	$ itap)
2980
2981	CALL VTe(VTinca)
2982	CALL VTb(VTphysiq)
2983	#endif
2984	END IF !config_inca /= 'none'
2985	c
2986	c Calculer les parametres optiques des nuages et quelques
2987	c parametres pour diagnostiques:
2988	c
2989
2990	IF (aerosol_couple) THEN
2991	mass_solu_aero(:,:) = ccm(:,:,1)
2992	mass_solu_aero_pi(:,:) = ccm(:,:,2)
2993	END IF
2994
2995	if (ok_newmicro) then
2996	CALL newmicro (paprs, pplay,ok_newmicro,
2997	. t_seri, cldliq, cldfra, cldtau, cldemi,
2998	. cldh, cldl, cldm, cldt, cldq,
2999	. flwp, fiwp, flwc, fiwc,
3000	e ok_aie,
3001	e mass_solu_aero, mass_solu_aero_pi,
3002	e bl95_b0, bl95_b1,
3003	s cldtaupi, re, fl, ref_liq, ref_ice)
3004	else
3005	CALL nuage (paprs, pplay,
3006	. t_seri, cldliq, cldfra, cldtau, cldemi,
3007	. cldh, cldl, cldm, cldt, cldq,
3008	e ok_aie,
3009	e mass_solu_aero, mass_solu_aero_pi,
3010	e bl95_b0, bl95_b1,
3011	s cldtaupi, re, fl)
3012
3013	endif
3014	c
3015	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3016	c
3017	IF (MOD(itaprad,radpas).EQ.0) THEN
3018
3019	DO i = 1, klon
3020	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
3021	. + falb1(i,is_lic) * pctsrf(i,is_lic)
3022	. + falb1(i,is_ter) * pctsrf(i,is_ter)
3023	. + falb1(i,is_sic) * pctsrf(i,is_sic)
3024	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
3025	. + falb2(i,is_lic) * pctsrf(i,is_lic)
3026	. + falb2(i,is_ter) * pctsrf(i,is_ter)
3027	. + falb2(i,is_sic) * pctsrf(i,is_sic)
3028	ENDDO
3029
3030	if (mydebug) then
3031	call writefield_phy('u_seri',u_seri,llm)
3032	call writefield_phy('v_seri',v_seri,llm)
3033	call writefield_phy('t_seri',t_seri,llm)
3034	call writefield_phy('q_seri',q_seri,llm)
3035	endif
3036
3037	IF (aerosol_couple) THEN
3038	#ifdef INCA
3039	CALL radlwsw_inca
3040	e (kdlon,kflev,dist, rmu0, fract, solaire,
3041	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
3042	e wo(:, :, 1),
3043	e cldfra, cldemi, cldtau,
3044	s heat,heat0,cool,cool0,radsol,albpla,
3045	s topsw,toplw,solsw,sollw,
3046	s sollwdown,
3047	s topsw0,toplw0,solsw0,sollw0,
3048	s lwdn0, lwdn, lwup0, lwup,
3049	s swdn0, swdn, swup0, swup,
3050	e ok_ade, ok_aie,
3051	e tau_aero, piz_aero, cg_aero,
3052	s topswad_aero, solswad_aero,
3053	s topswad0_aero, solswad0_aero,
3054	s topsw_aero, topsw0_aero,
3055	s solsw_aero, solsw0_aero,
3056	e cldtaupi,
3057	s topswai_aero, solswai_aero)
3058
3059	#endif
3060	ELSE
3061
3062	CALL radlwsw
3063	e (dist, rmu0, fract,
3064	e paprs, pplay,zxtsol,albsol1, albsol2,
3065	e t_seri,q_seri,wo,
3066	e cldfra, cldemi, cldtau,
3067	e ok_ade, ok_aie,
3068	e tau_aero, piz_aero, cg_aero,
3069	e cldtaupi,new_aod,
3070	e zqsat, flwc, fiwc,
3071	s heat,heat0,cool,cool0,radsol,albpla,
3072	s topsw,toplw,solsw,sollw,
3073	s sollwdown,
3074	s topsw0,toplw0,solsw0,sollw0,
3075	s lwdn0, lwdn, lwup0, lwup,
3076	s swdn0, swdn, swup0, swup,
3077	s topswad_aero, solswad_aero,
3078	s topswai_aero, solswai_aero,
3079	o topswad0_aero, solswad0_aero,
3080	o topsw_aero, topsw0_aero,
3081	o solsw_aero, solsw0_aero,
3082	o topswcf_aero, solswcf_aero)
3083
3084
3085	ENDIF ! aerosol_couple
3086	itaprad = 0
3087	ENDIF ! MOD(itaprad,radpas)
3088	itaprad = itaprad + 1
3089
3090	IF (iflag_radia.eq.0) THEN
3091	IF (prt_level.ge.9) THEN
3092	PRINT *,'--------------------------------------------------'
3093	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3094	PRINT *,'>>>> heat et cool mis a zero '
3095	PRINT *,'--------------------------------------------------'
3096	END IF
3097	heat=0.
3098	cool=0.
3099	END IF
3100
3101	c
3102	c Ajouter la tendance des rayonnements (tous les pas)
3103	c
3104	DO k = 1, klev
3105	DO i = 1, klon
3106	t_seri(i,k) = t_seri(i,k)
3107	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
3108	ENDDO
3109	ENDDO
3110	c
3111	if (mydebug) then
3112	call writefield_phy('u_seri',u_seri,llm)
3113	call writefield_phy('v_seri',v_seri,llm)
3114	call writefield_phy('t_seri',t_seri,llm)
3115	call writefield_phy('q_seri',q_seri,llm)
3116	endif
3117
3118	cIM
3119	IF (ip_ebil_phy.ge.2) THEN
3120	ztit='after rad'
3121	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3122	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3123	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3124	call diagphy(airephy,ztit,ip_ebil_phy
3125	e , topsw, toplw, solsw, sollw, zero_v
3126	e , zero_v, zero_v, zero_v, ztsol
3127	e , d_h_vcol, d_qt, d_ec
3128	s , fs_bound, fq_bound )
3129	END IF
3130	c
3131	c
3132	c Calculer l'hydrologie de la surface
3133	c
3134	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3135	c . agesno, ftsol,fqsurf,fsnow, ruis)
3136	c
3137
3138	c
3139	c Calculer le bilan du sol et la derive de temperature (couplage)
3140	c
3141	DO i = 1, klon
3142	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3143	c a la demande de JLD
3144	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3145	ENDDO
3146	c
3147	cmoddeblott(jan95)
3148	c Appeler le programme de parametrisation de l'orographie
3149	c a l'echelle sous-maille:
3150	c
3151	IF (ok_orodr) THEN
3152	c
3153	c selection des points pour lesquels le shema est actif:
3154	igwd=0
3155	DO i=1,klon
3156	itest(i)=0
3157	c IF ((zstd(i).gt.10.0)) THEN
3158	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3159	itest(i)=1
3160	igwd=igwd+1
3161	idx(igwd)=i
3162	ENDIF
3163	ENDDO
3164	c igwdim=MAX(1,igwd)
3165	c
3166	IF (ok_strato) THEN
3167
3168	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3169	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3170	e igwd,idx,itest,
3171	e t_seri, u_seri, v_seri,
3172	s zulow, zvlow, zustrdr, zvstrdr,
3173	s d_t_oro, d_u_oro, d_v_oro)
3174
3175	ELSE
3176	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3177	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3178	e igwd,idx,itest,
3179	e t_seri, u_seri, v_seri,
3180	s zulow, zvlow, zustrdr, zvstrdr,
3181	s d_t_oro, d_u_oro, d_v_oro)
3182	ENDIF
3183	c
3184	c ajout des tendances
3185	!-----------------------------------------------------------------------------------------
3186	! ajout des tendances de la trainee de l'orographie
3187	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3188	!-----------------------------------------------------------------------------------------
3189	c
3190	ENDIF ! fin de test sur ok_orodr
3191	c
3192	if (mydebug) then
3193	call writefield_phy('u_seri',u_seri,llm)
3194	call writefield_phy('v_seri',v_seri,llm)
3195	call writefield_phy('t_seri',t_seri,llm)
3196	call writefield_phy('q_seri',q_seri,llm)
3197	endif
3198
3199	IF (ok_orolf) THEN
3200	c
3201	c selection des points pour lesquels le shema est actif:
3202	igwd=0
3203	DO i=1,klon
3204	itest(i)=0
3205	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3206	itest(i)=1
3207	igwd=igwd+1
3208	idx(igwd)=i
3209	ENDIF
3210	ENDDO
3211	c igwdim=MAX(1,igwd)
3212	c
3213	IF (ok_strato) THEN
3214
3215	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3216	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3217	e igwd,idx,itest,
3218	e t_seri, u_seri, v_seri,
3219	s zulow, zvlow, zustrli, zvstrli,
3220	s d_t_lif, d_u_lif, d_v_lif )
3221
3222	ELSE
3223	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3224	e rlat,zmea,zstd,zpic,
3225	e itest,
3226	e t_seri, u_seri, v_seri,
3227	s zulow, zvlow, zustrli, zvstrli,
3228	s d_t_lif, d_u_lif, d_v_lif)
3229	ENDIF
3230	c
3231	!-----------------------------------------------------------------------------------------
3232	! ajout des tendances de la portance de l'orographie
3233	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3234	!-----------------------------------------------------------------------------------------
3235	c
3236	ENDIF ! fin de test sur ok_orolf
3237	C HINES GWD PARAMETRIZATION
3238
3239	IF (ok_hines) then
3240
3241	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3242	i rlat,t_seri,u_seri,v_seri,
3243	o zustrhi,zvstrhi,
3244	o d_t_hin, d_u_hin, d_v_hin)
3245	c
3246	c ajout des tendances
3247	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3248
3249	ENDIF
3250	c
3251
3252	c
3253	cIM cf. FLott BEG
3254	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3255
3256	if (mydebug) then
3257	call writefield_phy('u_seri',u_seri,llm)
3258	call writefield_phy('v_seri',v_seri,llm)
3259	call writefield_phy('t_seri',t_seri,llm)
3260	call writefield_phy('q_seri',q_seri,llm)
3261	endif
3262
3263	DO i = 1, klon
3264	zustrph(i)=0.
3265	zvstrph(i)=0.
3266	ENDDO
3267	DO k = 1, klev
3268	DO i = 1, klon
3269	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3270	c (paprs(i,k)-paprs(i,k+1))/rg
3271	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3272	c (paprs(i,k)-paprs(i,k+1))/rg
3273	ENDDO
3274	ENDDO
3275	c
3276	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3277	c
3278	IF (is_sequential) THEN
3279
3280	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3281	C ra,rg,romega,
3282	C rlat,rlon,pphis,
3283	C zustrdr,zustrli,zustrph,
3284	C zvstrdr,zvstrli,zvstrph,
3285	C paprs,u,v,
3286	C aam, torsfc)
3287	ENDIF
3288	cIM cf. FLott END
3289	cIM
3290	IF (ip_ebil_phy.ge.2) THEN
3291	ztit='after orography'
3292	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3293	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3294	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3295	call diagphy(airephy,ztit,ip_ebil_phy
3296	e , zero_v, zero_v, zero_v, zero_v, zero_v
3297	e , zero_v, zero_v, zero_v, ztsol
3298	e , d_h_vcol, d_qt, d_ec
3299	s , fs_bound, fq_bound )
3300	END IF
3301	c
3302	c
3303	!====================================================================
3304	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3305	!====================================================================
3306	! Abderrahmane 24.08.09
3307
3308	IF (ok_cosp) THEN
3309	! adeclarer
3310	#ifdef CPP_COSP
3311	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3312
3313	print*,'freq_cosp',freq_cosp
3314	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3315	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3316	! s ref_liq,ref_ice
3317	call phys_cosp(itap,dtime,freq_cosp,
3318	$ ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
3319	$ ecrit_mth,ecrit_day,ecrit_hf,
3320	$ klon,klev,rlon,rlat,presnivs,overlap,
3321	$ ref_liq,ref_ice,
3322	$ pctsrf(:,is_ter)+pctsrf(:,is_lic),
3323	$ zu10m,zv10m,
3324	$ zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3325	$ qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3326	$ prfl(:,1:klev),psfl(:,1:klev),
3327	$ pmflxr(:,1:klev),pmflxs(:,1:klev),
3328	$ mr_ozone,cldtau, cldemi)
3329
3330	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3331	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3332	! M clMISR,
3333	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3334	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3335
3336	ENDIF
3337
3338	#endif
3339	ENDIF !ok_cosp
3340	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3341	cAA
3342	cAA Installation de l'interface online-offline pour traceurs
3343	cAA
3344	c====================================================================
3345	c Calcul des tendances traceurs
3346	c====================================================================
3347	C
3348
3349	call phytrac (
3350	I itap, days_elapsed+1, jH_cur, debut,
3351	I lafin, dtime, u, v, t,
3352	I paprs, pplay, pmfu, pmfd,
3353	I pen_u, pde_u, pen_d, pde_d,
3354	I cdragh, coefh, fm_therm, entr_therm,
3355	I u1, v1, ftsol, pctsrf,
3356	I rlat, frac_impa, frac_nucl,rlon,
3357	I presnivs, pphis, pphi, albsol1,
3358	I qx(:,:,ivap),rhcl, cldfra, rneb,
3359	I diafra, cldliq, itop_con, ibas_con,
3360	I pmflxr, pmflxs, prfl, psfl,
3361	I da, phi, mp, upwd,
3362	I dnwd, aerosol_couple, flxmass_w,
3363	I tau_aero, piz_aero, cg_aero, ccm,
3364	I rfname,
3365	O tr_seri)
3366
3367	IF (offline) THEN
3368
3369	IF (prt_level.ge.9)
3370	$ print*,'Attention on met a 0 les thermiques pour phystoke'
3371	call phystokenc (
3372	I nlon,klev,pdtphys,rlon,rlat,
3373	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3374	I fm_therm,entr_therm,
3375	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3376	I frac_impa, frac_nucl,
3377	I pphis,airephy,dtime,itap,
3378	I rlon,rlat,qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3379
3380
3381	ENDIF
3382
3383	c
3384	c Calculer le transport de l'eau et de l'energie (diagnostique)
3385	c
3386	CALL transp (paprs,zxtsol,
3387	e t_seri, q_seri, u_seri, v_seri, zphi,
3388	s ve, vq, ue, uq)
3389	c
3390	cIM global posePB BEG
3391	IF(1.EQ.0) THEN
3392	c
3393	CALL transp_lay (paprs,zxtsol,
3394	e t_seri, q_seri, u_seri, v_seri, zphi,
3395	s ve_lay, vq_lay, ue_lay, uq_lay)
3396	c
3397	ENDIF !(1.EQ.0) THEN
3398	cIM global posePB END
3399	c Accumuler les variables a stocker dans les fichiers histoire:
3400	c
3401	c+jld ec_conser
3402	DO k = 1, klev
3403	DO i = 1, klon
3404	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3405	d_t_ec(i,k)=0.5/ZRCPD
3406	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3407	ENDDO
3408	ENDDO
3409
3410	DO k = 1, klev
3411	DO i = 1, klon
3412	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3413	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3414	END DO
3415	END DO
3416	c-jld ec_conser
3417	cIM
3418	IF (ip_ebil_phy.ge.1) THEN
3419	ztit='after physic'
3420	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3421	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3422	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3423	C Comme les tendances de la physique sont ajoute dans la dynamique,
3424	C on devrait avoir que la variation d'entalpie par la dynamique
3425	C est egale a la variation de la physique au pas de temps precedent.
3426	C Donc la somme de ces 2 variations devrait etre nulle.
3427
3428	call diagphy(airephy,ztit,ip_ebil_phy
3429	e , topsw, toplw, solsw, sollw, sens
3430	e , evap, rain_fall, snow_fall, ztsol
3431	e , d_h_vcol, d_qt, d_ec
3432	s , fs_bound, fq_bound )
3433	C
3434	d_h_vcol_phy=d_h_vcol
3435	C
3436	END IF
3437	C
3438	c=======================================================================
3439	c SORTIES
3440	c=======================================================================
3441
3442	cIM Interpolation sur les niveaux de pression du NMC
3443	c -------------------------------------------------
3444	c
3445	#include "calcul_STDlev.h"
3446	twriteSTD(:,:,1)=tsumSTD(:,:,1)
3447	qwriteSTD(:,:,1)=qsumSTD(:,:,1)
3448	rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
3449	phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
3450	uwriteSTD(:,:,1)=usumSTD(:,:,1)
3451	vwriteSTD(:,:,1)=vsumSTD(:,:,1)
3452	wwriteSTD(:,:,1)=wsumSTD(:,:,1)
3453
3454	twriteSTD(:,:,2)=tsumSTD(:,:,2)
3455	qwriteSTD(:,:,2)=qsumSTD(:,:,2)
3456	rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
3457	phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
3458	uwriteSTD(:,:,2)=usumSTD(:,:,2)
3459	vwriteSTD(:,:,2)=vsumSTD(:,:,2)
3460	wwriteSTD(:,:,2)=wsumSTD(:,:,2)
3461
3462	twriteSTD(:,:,3)=tlevSTD(:,:)
3463	qwriteSTD(:,:,3)=qlevSTD(:,:)
3464	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3465	phiwriteSTD(:,:,3)=philevSTD(:,:)
3466	uwriteSTD(:,:,3)=ulevSTD(:,:)
3467	vwriteSTD(:,:,3)=vlevSTD(:,:)
3468	wwriteSTD(:,:,3)=wlevSTD(:,:)
3469
3470	twriteSTD(:,:,4)=tlevSTD(:,:)
3471	qwriteSTD(:,:,4)=qlevSTD(:,:)
3472	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3473	phiwriteSTD(:,:,4)=philevSTD(:,:)
3474	uwriteSTD(:,:,4)=ulevSTD(:,:)
3475	vwriteSTD(:,:,4)=vlevSTD(:,:)
3476	wwriteSTD(:,:,4)=wlevSTD(:,:)
3477	c
3478	cIM initialisation 5eme fichier de sortie
3479	cIM ajoute 5eme niveau 170310 BEG
3480	twriteSTD(:,:,5)=tlevSTD(:,:)
3481	qwriteSTD(:,:,5)=qlevSTD(:,:)
3482	rhwriteSTD(:,:,5)=rhlevSTD(:,:)
3483	phiwriteSTD(:,:,5)=philevSTD(:,:)
3484	uwriteSTD(:,:,5)=ulevSTD(:,:)
3485	vwriteSTD(:,:,5)=vlevSTD(:,:)
3486	wwriteSTD(:,:,5)=wlevSTD(:,:)
3487	cIM for NMC files
3488	DO n=1, nlevSTD3
3489	DO k=1, nlevSTD
3490	if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
3491	twriteSTD3(:,n)=tlevSTD(:,k)
3492	qwriteSTD3(:,n)=qlevSTD(:,k)
3493	rhwriteSTD3(:,n)=rhlevSTD(:,k)
3494	phiwriteSTD3(:,n)=philevSTD(:,k)
3495	uwriteSTD3(:,n)=ulevSTD(:,k)
3496	vwriteSTD3(:,n)=vlevSTD(:,k)
3497	wwriteSTD3(:,n)=wlevSTD(:,k)
3498	endif !rlevSTD3(n).EQ.rlevSTD(k)
3499	ENDDO
3500	ENDDO
3501	c
3502	DO n=1, nlevSTD8
3503	DO k=1, nlevSTD
3504	if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
3505	tnondefSTD8(:,n)=tnondef(:,k,2)
3506	twriteSTD8(:,n)=tsumSTD(:,k,2)
3507	qwriteSTD8(:,n)=qsumSTD(:,k,2)
3508	rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
3509	phiwriteSTD8(:,n)=phisumSTD(:,k,2)
3510	uwriteSTD8(:,n)=usumSTD(:,k,2)
3511	vwriteSTD8(:,n)=vsumSTD(:,k,2)
3512	wwriteSTD8(:,n)=wsumSTD(:,k,2)
3513	endif !rlevSTD8(n).EQ.rlevSTD(k)
3514	ENDDO
3515	ENDDO
3516	c
3517	c slp sea level pressure
3518	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3519	c
3520	ccc prw = eau precipitable
3521	DO i = 1, klon
3522	prw(i) = 0.
3523	DO k = 1, klev
3524	prw(i) = prw(i) +
3525	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3526	ENDDO
3527	ENDDO
3528	c
3529	cIM initialisation + calculs divers diag AMIP2
3530	c
3531	#include "calcul_divers.h"
3532	c
3533	IF (config_inca /= 'none') THEN
3534	#ifdef INCA
3535	CALL VTe(VTphysiq)
3536	CALL VTb(VTinca)
3537
3538	CALL chemhook_end (
3539	$ dtime,
3540	$ pplay,
3541	$ t_seri,
3542	$ tr_seri,
3543	$ nbtr,
3544	$ paprs,
3545	$ q_seri,
3546	$ airephy,
3547	$ pphi,
3548	$ pphis,
3549	$ zx_rh)
3550
3551	CALL VTe(VTinca)
3552	CALL VTb(VTphysiq)
3553	#endif
3554	END IF
3555
3556	c=============================================================
3557	c
3558	c Convertir les incrementations en tendances
3559	c
3560	if (mydebug) then
3561	call writefield_phy('u_seri',u_seri,llm)
3562	call writefield_phy('v_seri',v_seri,llm)
3563	call writefield_phy('t_seri',t_seri,llm)
3564	call writefield_phy('q_seri',q_seri,llm)
3565	endif
3566
3567	DO k = 1, klev
3568	DO i = 1, klon
3569	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3570	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3571	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3572	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3573	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3574	ENDDO
3575	ENDDO
3576	c
3577	IF (nqtot.GE.3) THEN
3578	DO iq = 3, nqtot
3579	DO k = 1, klev
3580	DO i = 1, klon
3581	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3582	ENDDO
3583	ENDDO
3584	ENDDO
3585	ENDIF
3586	c
3587	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3588	cIM global posePB#include "write_bilKP_ins.h"
3589	cIM global posePB#include "write_bilKP_ave.h"
3590	c
3591
3592	c Sauvegarder les valeurs de t et q a la fin de la physique:
3593	c
3594	DO k = 1, klev
3595	DO i = 1, klon
3596	u_ancien(i,k) = u_seri(i,k)
3597	v_ancien(i,k) = v_seri(i,k)
3598	t_ancien(i,k) = t_seri(i,k)
3599	q_ancien(i,k) = q_seri(i,k)
3600	ENDDO
3601	ENDDO
3602	c
3603	!==========================================================================
3604	! Sorties des tendances pour un point particulier
3605	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3606	! pour le debug
3607	! La valeur de igout est attribuee plus haut dans le programme
3608	!==========================================================================
3609
3610	if (prt_level.ge.1) then
3611	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3612	write(lunout,*)
3613	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3614	write(lunout,*)
3615	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3616	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3617	s pctsrf(igout,is_sic)
3618	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3619	do k=1,klev
3620	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3621	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3622	s d_t_eva(igout,k)
3623	enddo
3624	write(lunout,*) 'cool,heat'
3625	do k=1,klev
3626	write(lunout,*) cool(igout,k),heat(igout,k)
3627	enddo
3628
3629	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3630	do k=1,klev
3631	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3632	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3633	enddo
3634
3635	write(lunout,*) 'd_ps ',d_ps(igout)
3636	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3637	do k=1,klev
3638	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3639	s d_qx(igout,k,1),d_qx(igout,k,2)
3640	enddo
3641	endif
3642
3643	!==========================================================================
3644
3645	c============================================================
3646	c Calcul de la temperature potentielle
3647	c============================================================
3648	DO k = 1, klev
3649	DO i = 1, klon
3650	cJYG/IM theta en debut du pas de temps
3651	cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3652	cJYG/IM theta en fin de pas de temps de physique
3653	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
3654	ENDDO
3655	ENDDO
3656	c
3657
3658	c 22.03.04 BEG
3659	c=============================================================
3660	c Ecriture des sorties
3661	c=============================================================
3662	#ifdef CPP_IOIPSL
3663
3664	c Recupere des varibles calcule dans differents modules
3665	c pour ecriture dans histxxx.nc
3666
3667	! Get some variables from module fonte_neige_mod
3668	CALL fonte_neige_get_vars(pctsrf,
3669	. zxfqcalving, zxfqfonte, zxffonte)
3670
3671
3672	#include "phys_output_write.h"
3673
3674	#ifdef histISCCP
3675	#include "write_histISCCP.h"
3676	#endif
3677
3678	#ifdef histNMC
3679	#include "write_histhfNMC.h"
3680	#include "write_histdayNMC.h"
3681	#include "write_histmthNMC.h"
3682	#endif
3683
3684	#include "write_histday_seri.h"
3685
3686	#include "write_paramLMDZ_phy.h"
3687
3688	#endif
3689
3690	c 22.03.04 END
3691	c
3692	c====================================================================
3693	c Si c'est la fin, il faut conserver l'etat de redemarrage
3694	c====================================================================
3695	c
3696
3697
3698	IF (lafin) THEN
3699	itau_phy = itau_phy + itap
3700	CALL phyredem ("restartphy.nc")
3701	! open(97,form="unformatted",file="finbin")
3702	! write(97) u_seri,v_seri,t_seri,q_seri
3703	! close(97)
3704	C$OMP MASTER
3705	if (read_climoz >= 1) then
3706	if (is_mpi_root) then
3707	call nf95_close(ncid_climoz)
3708	end if
3709	deallocate(press_climoz) ! pointer
3710	end if
3711	C$OMP END MASTER
3712	ENDIF
3713
3714	! first=.false.
3715
3716	RETURN
3717	END
3718	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3719	IMPLICIT none
3720	c
3721	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3722	c la conservation de l'eau
3723	c
3724	#include "YOMCST.h"
3725	INTEGER klon,klev
3726	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3727	REAL aire(klon)
3728	REAL qtotal, zx, qcheck
3729	INTEGER i, k
3730	c
3731	zx = 0.0
3732	DO i = 1, klon
3733	zx = zx + aire(i)
3734	ENDDO
3735	qtotal = 0.0
3736	DO k = 1, klev
3737	DO i = 1, klon
3738	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3739	. *(paprs(i,k)-paprs(i,k+1))/RG
3740	ENDDO
3741	ENDDO
3742	c
3743	qcheck = qtotal/zx
3744	c
3745	RETURN
3746	END
3747	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3748	IMPLICIT none
3749	c
3750	c Tranformer une variable de la grille physique a
3751	c la grille d'ecriture
3752	c
3753	INTEGER nfield,nlon,iim,jjmp1, jjm
3754	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3755	c
3756	INTEGER i, n, ig
3757	c
3758	jjm = jjmp1 - 1
3759	DO n = 1, nfield
3760	DO i=1,iim
3761	ecrit(i,n) = fi(1,n)
3762	ecrit(i+jjm*iim,n) = fi(nlon,n)
3763	ENDDO
3764	DO ig = 1, nlon - 2
3765	ecrit(iim+ig,n) = fi(1+ig,n)
3766	ENDDO
3767	ENDDO
3768	RETURN
3769	END

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

Download in other formats:

Original Format