Context Navigation

physiq.F @ 1536

Last change on this file since 1536 was 1536, checked in by musat, 13 years ago
Petites corrections pour le cas ou l'on veut pas sortir de fichier histstn (et qu'on a pas les fichiers initiaux correspondants). IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 127.0 KB

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

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

Download in other formats:

Original Format