Context Navigation

physiq.F @ 1226

Last change on this file since 1226 was 1226, checked in by jghattas, 15 years ago

La fonction ioget_year_len n'existe pas dans les tags officiels
d'IOIPSL, uniquement dans le trunk. LMDZ n'compilant plus, j'ai
commente l'appel a cette fonction et temporairement desactive l'option
read_climoz qui l'utilise.

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

File size: 115.1 KB

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

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

Download in other formats:

Original Format