Context Navigation

physiq.F @ 1869

Last change on this file since 1869 was 1717, checked in by Ehouarn Millour, 12 years ago
Added "arch" files for Ada (using dynamic libraries for NetCDF, you must have export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/smplocal/pub/NetCDF/4.1.3/lib:/smplocal/pub/HDF5/1.8.9/seq/lib in your .bashrc or .bash_login or in your job to run). Also updated some sources so that gcm bench runs in "debug" mode (note that all these changes are minor and have already been implemented in LMDZ5 trunk). EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 127.2 KB

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

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

Download in other formats:

Original Format