Context Navigation

physiq.F @ 1086

Last change on this file since 1086 was 740, checked in by lmdzadmin, 18 years ago
Correction bogues: les ecrit_ sont des REALs lus dans conf_phys.F90 en nombre de jours sauf pour ecrit_ins et ecrit_tra en secondes! Les ecrit_ sont initialises dans conf_phys.F90 et peuvent etre modifies dans physiq.def. IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.2 KB

Line
1	!
2	! $Header$
3	!
4	c
5	SUBROUTINE physiq (nlon,nlev,nqmax,
6	. debut,lafin,rjourvrai,gmtime,pdtphys,
7	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
8	. u,v,t,qx,
9	. omega,
10	#ifdef INCA
11	. flxmass_w,
12	#endif
13	. d_u, d_v, d_t, d_qx, d_ps
14	cIM Amip2
15	. , dudyn
16	. , PVteta)
17
18	USE ioipsl
19	USE histcom
20	#ifdef INCA
21	USE chemshut
22	USE species_names
23	#ifdef INCA_CH4
24	! USE obs_pos
25	#endif
26	#endif
27	IMPLICIT none
28	c======================================================================
29	c
30	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
31	c
32	c Objet: Moniteur general de la physique du modele
33	cAA Modifications quant aux traceurs :
34	cAA - uniformisation des parametrisations ds phytrac
35	cAA - stockage des moyennes des champs necessaires
36	cAA en mode traceur off-line
37	c======================================================================
38	c CLEFS CPP POUR LES IO
39	c =====================
40	#define histins
41	#define histhf
42	#define histday
43	#define histmth
44	#define histISCCP
45	#define histmthNMC
46	c======================================================================
47	c modif ( P. Le Van , 12/10/98 )
48	c
49	c Arguments:
50	c
51	c nlon----input-I-nombre de points horizontaux
52	c nlev----input-I-nombre de couches verticales
53	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
54	c debut---input-L-variable logique indiquant le premier passage
55	c lafin---input-L-variable logique indiquant le dernier passage
56	c rjour---input-R-numero du jour de l'experience
57	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
58	c pdtphys-input-R-pas d'integration pour la physique (seconde)
59	c paprs---input-R-pression pour chaque inter-couche (en Pa)
60	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
61	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
62	c pphis---input-R-geopotentiel du sol
63	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
64	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
65	c v-------input-R-vitesse Y (de S a N) en m/s
66	c t-------input-R-temperature (K)
67	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
68	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
69	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
70	c omega---input-R-vitesse verticale en Pa/s
71	c
72	c d_u-----output-R-tendance physique de "u" (m/s/s)
73	c d_v-----output-R-tendance physique de "v" (m/s/s)
74	c d_t-----output-R-tendance physique de "t" (K/s)
75	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
76	c d_ps----output-R-tendance physique de la pression au sol
77	cIM
78	c PVteta--output-R-vorticite potentielle a des thetas constantes
79	c======================================================================
80	#include "dimensions.h"
81	integer jjmp1
82	parameter (jjmp1=jjm+1-1/jjm)
83	#include "dimphy.h"
84	#include "regdim.h"
85	#include "indicesol.h"
86	#include "dimsoil.h"
87	#include "clesphys.h"
88	#include "control.h"
89	#include "logic.h"
90	#include "temps.h"
91	#include "comgeomphy.h"
92	#include "advtrac.h"
93	#include "iniprint.h"
94	#include "thermcell.h"
95	c======================================================================
96	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
97	PARAMETER (ok_cvl=.TRUE.)
98	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
99	PARAMETER (ok_gust=.FALSE.)
100	c======================================================================
101	LOGICAL check ! Verifier la conservation du modele en eau
102	PARAMETER (check=.FALSE.)
103	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
104	PARAMETER (ok_stratus=.FALSE.)
105	c======================================================================
106	c Parametres lies au coupleur OASIS:
107	#include "oasis.h"
108	INTEGER,SAVE :: npas, nexca
109	logical rnpb
110	#ifdef INCA
111	parameter(rnpb=.false.)
112	#else
113	parameter(rnpb=.true.)
114	#endif
115	c ocean = type de modele ocean a utiliser: force, slab, couple
116	character*6 ocean
117	SAVE ocean
118
119	c parameter (ocean = 'force ')
120	c parameter (ocean = 'couple')
121	logical ok_ocean
122	SAVE ok_ocean
123	c
124	cIM "slab" ocean
125	REAL tslab(klon) !Temperature du slab-ocean
126	REAL seaice(klon) !glace de mer (kg/m2)
127	REAL fluxo(klon) !flux turbulents ocean-glace de mer
128	REAL fluxg(klon) !flux turbulents ocean-atmosphere
129	c
130	c======================================================================
131	c Clef controlant l'activation du cycle diurne:
132	ccc LOGICAL cycle_diurne
133	ccc PARAMETER (cycle_diurne=.FALSE.)
134	c======================================================================
135	c Modele thermique du sol, a activer pour le cycle diurne:
136	ccc LOGICAL soil_model
137	ccc PARAMETER (soil_model=.FALSE.)
138	logical ok_veget
139	save ok_veget
140	c parameter (ok_veget = .true.)
141	c parameter (ok_veget = .false.)
142	c======================================================================
143	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
144	c le calcul du rayonnement est celle apres la precipitation des nuages.
145	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
146	c la condensation et la precipitation. Cette cle augmente les impacts
147	c radiatifs des nuages.
148	ccc LOGICAL new_oliq
149	ccc PARAMETER (new_oliq=.FALSE.)
150	c======================================================================
151	c Clefs controlant deux parametrisations de l'orographie:
152	cc LOGICAL ok_orodr
153	ccc PARAMETER (ok_orodr=.FALSE.)
154	ccc LOGICAL ok_orolf
155	ccc PARAMETER (ok_orolf=.FALSE.)
156	c======================================================================
157	LOGICAL ok_journe ! sortir le fichier journalier
158	save ok_journe
159	c PARAMETER (ok_journe=.true.)
160	c
161	LOGICAL ok_mensuel ! sortir le fichier mensuel
162	save ok_mensuel
163	c PARAMETER (ok_mensuel=.true.)
164	c
165	LOGICAL ok_instan ! sortir le fichier instantane
166	save ok_instan
167	c PARAMETER (ok_instan=.true.)
168	c
169	LOGICAL ok_region ! sortir le fichier regional
170	PARAMETER (ok_region=.FALSE.)
171	c======================================================================
172	c pour phsystoke avec thermiques
173	REAL fm_therm(klon,klev+1)
174	REAL entr_therm(klon,klev)
175	real q2(klon,klev+1,nbsrf)
176	save q2
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	INTEGER nqmax
191	REAL rjourvrai
192	REAL gmtime
193	REAL pdtphys
194	LOGICAL debut, lafin
195	REAL paprs(klon,klev+1)
196	REAL pplay(klon,klev)
197	REAL pphi(klon,klev)
198	REAL pphis(klon)
199	REAL presnivs(klev)
200	REAL znivsig(klev)
201	REAL zsurf(nbsrf)
202	cIM
203	INTEGER kinv
204	real pir
205	cMI
206	REAL u(klon,klev)
207	REAL v(klon,klev)
208	REAL t(klon,klev)
209	REAL qx(klon,klev,nqmax)
210
211	REAL t_ancien(klon,klev), q_ancien(klon,klev)
212	SAVE t_ancien, q_ancien
213	LOGICAL ancien_ok
214	SAVE ancien_ok
215
216	REAL d_t_dyn(klon,klev)
217	REAL d_q_dyn(klon,klev)
218
219	REAL omega(klon,klev)
220
221	#ifdef INCA
222	REAL flxmass_w(klon,klev)
223	#endif
224	REAL d_u(klon,klev)
225	REAL d_v(klon,klev)
226	REAL d_t(klon,klev)
227	REAL d_qx(klon,klev,nqmax)
228	REAL d_ps(klon)
229	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
230	c
231	cIM Amip2 PV a theta constante
232	c
233	INTEGER nbteta
234	PARAMETER(nbteta=3)
235	CHARACTER*3 ctetaSTD(nbteta)
236	DATA ctetaSTD/'350','380','405'/
237	REAL rtetaSTD(nbteta)
238	DATA rtetaSTD/350., 380., 405./
239	c
240	REAL PVteta(klon,nbteta)
241	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
242	c
243	cMI Amip2 PV a theta constante
244
245	INTEGER klevp1, klevm1
246	PARAMETER(klevp1=klev+1,klevm1=klev-1)
247	#include "raddim.h"
248	c
249	cIM 080304 REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
250	REAL swdn0(klon,klevp1), swdn(klon,klevp1)
251	REAL swup0(klon,klevp1), swup(klon,klevp1)
252	SAVE swdn0 , swdn, swup0, swup
253	c
254	REAL SWdn200clr(klon), SWdn200(klon)
255	REAL SWup200clr(klon), SWup200(klon)
256	SAVE SWdn200clr, SWdn200, SWup200clr, SWup200
257	c
258	REAL lwdn0(klon,klevp1), lwdn(klon,klevp1)
259	REAL lwup0(klon,klevp1), lwup(klon,klevp1)
260	SAVE lwdn0 , lwdn, lwup0, lwup
261	c
262	REAL LWdn200clr(klon), LWdn200(klon)
263	REAL LWup200clr(klon), LWup200(klon)
264	SAVE LWdn200clr, LWdn200, LWup200clr, LWup200
265	c
266	REAL LWdnTOA(klon), LWdnTOAclr(klon)
267	SAVE LWdnTOA, LWdnTOAclr
268	c
269	cIM Amip2
270	c variables a une pression donnee
271	c
272	integer nlevSTD
273	PARAMETER(nlevSTD=17)
274	real rlevSTD(nlevSTD)
275	DATA rlevSTD/100000., 92500., 85000., 70000.,
276	.60000., 50000., 40000., 30000., 25000., 20000.,
277	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
278	CHARACTER*4 clevSTD(nlevSTD)
279	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
280	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
281	.'70 ','50 ','30 ','20 ','10 '/
282	c
283	CHARACTER*3 bb2
284	CHARACTER*2 bb3
285	c
286	real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
287	real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
288	real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
289	real wlevSTD(klon,nlevSTD)
290	c
291	c nout : niveau de output des variables a une pression donnee
292	INTEGER nout
293	PARAMETER(nout=3) !nout=1 : day; =2 : mth; =3 : NMC
294	c
295	REAL tsumSTD(klon,nlevSTD,nout)
296	REAL usumSTD(klon,nlevSTD,nout), vsumSTD(klon,nlevSTD,nout)
297	REAL wsumSTD(klon,nlevSTD,nout), phisumSTD(klon,nlevSTD,nout)
298	REAL qsumSTD(klon,nlevSTD,nout), rhsumSTD(klon,nlevSTD,nout)
299	c
300	SAVE tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD,
301	. qsumSTD, rhsumSTD
302	c
303	logical oknondef(klon,nlevSTD,nout)
304	real tnondef(klon,nlevSTD,nout)
305	save tnondef
306	c
307	c les produits uvSTD, vqSTD, .., T2STD sont calcules
308	c a partir des valeurs instantannees toutes les 6 h
309	c qui sont moyennees sur le mois
310	c
311	real uvSTD(klon,nlevSTD)
312	real vqSTD(klon,nlevSTD)
313	real vTSTD(klon,nlevSTD)
314	real wqSTD(klon,nlevSTD)
315	c
316	real uvsumSTD(klon,nlevSTD,nout)
317	real vqsumSTD(klon,nlevSTD,nout)
318	real vTsumSTD(klon,nlevSTD,nout)
319	real wqsumSTD(klon,nlevSTD,nout)
320	c
321	real vphiSTD(klon,nlevSTD)
322	real wTSTD(klon,nlevSTD)
323	real u2STD(klon,nlevSTD)
324	real v2STD(klon,nlevSTD)
325	real T2STD(klon,nlevSTD)
326	c
327	real vphisumSTD(klon,nlevSTD,nout)
328	real wTsumSTD(klon,nlevSTD,nout)
329	real u2sumSTD(klon,nlevSTD,nout)
330	real v2sumSTD(klon,nlevSTD,nout)
331	real T2sumSTD(klon,nlevSTD,nout)
332	c
333	SAVE uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD
334	SAVE vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD
335	cMI Amip2
336	c
337	#include "radepsi.h"
338	#include "radopt.h"
339	c
340	c
341	c prw: precipitable water
342	real prw(klon)
343
344	REAL convliq(klon,klev) ! eau liquide nuageuse convective
345	REAL convfra(klon,klev) ! fraction nuageuse convective
346
347	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
348	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
349	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
350	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
351
352	INTEGER kp1
353	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
354	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
355	REAL flwp(klon), fiwp(klon)
356	REAL flwc(klon,klev), fiwc(klon,klev)
357	REAL flwp_c(klon), fiwp_c(klon)
358	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
359	REAL flwp_s(klon), fiwp_s(klon)
360	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
361
362	cIM ISCCP simulator v3.4
363	c dans clesphys.h top_height, overlap
364	cv3.4
365	INTEGER debug, debugcol
366	INTEGER npoints
367	PARAMETER(npoints=klon)
368	c
369	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
370	INTEGER nregISCtot
371	PARAMETER(nregISCtot=1)
372	c
373	c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
374	c y compris pour 1 point
375	c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
376	c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
377	INTEGER imin_debut, nbpti
378	INTEGER jmin_debut, nbptj
379	c
380	REAL nbsunlit(nregISCtot,klon) !nbsunlit : moyenne de sunlit
381	INTEGER ncol, seed(klon)
382
383	cIM ajout seed_re (REAL)et seed_old (integer)
384	REAL seed_re(klon), aa
385	INTEGER seed_old(klon)
386	SAVE seed_old
387
388
389
390	c ncol = nb. de sous-colonnes pour chaque maille du GCM
391	c PARAMETER(ncol=100)
392	c PARAMETER(ncol=625)
393	c PARAMETER(ncol=10)
394	PARAMETER(ncol=25)
395	REAL tautab(0:255)
396	INTEGER invtau(-20:45000)
397	REAL emsfc_lw
398	PARAMETER(emsfc_lw=0.99)
399	c REAL ran0 ! type for random number fuction
400	c
401	REAL cldtot(klon,klev)
402	c variables de haut en bas pour le simulateur ISCCP
403	REAL dtau_s(klon,klev) !tau nuages startiformes
404	REAL dtau_c(klon,klev) !tau nuages convectifs
405	REAL dem_s(klon,klev) !emissivite nuages startiformes
406	REAL dem_c(klon,klev) !emissivite nuages convectifs
407	c
408	c variables de haut en bas pour le simulateur ISCCP
409	REAL pfull(klon,klev)
410	REAL phalf(klon,klev+1)
411	REAL qv(klon,klev)
412	REAL cc(klon,klev)
413	REAL conv(klon,klev)
414	REAL dtau_sH2B(klon,klev)
415	REAL dtau_cH2B(klon,klev)
416	REAL at(klon,klev)
417	REAL dem_sH2B(klon,klev)
418	REAL dem_cH2B(klon,klev)
419
420	c output from ISCCP simulator
421	REAL fq_isccp(klon,7,7)
422	REAL totalcldarea(klon)
423	REAL meanptop(klon)
424	REAL meantaucld(klon)
425	REAL boxtau(klon,ncol)
426	REAL boxptop(klon,ncol)
427	c
428	INTEGER l, kmax, lmax
429	PARAMETER(kmax=8, lmax=8)
430	INTEGER kmaxm1, lmaxm1
431	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
432	INTEGER iimx7, jjmx7, jjmp1x7
433	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1,
434	.jjmp1x7=jjmp1*lmaxm1)
435	c
436	INTEGER iw, iwmax
437	REAL wmin, pas_w
438	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
439	PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
440	REAL o500(klon)
441	c
442
443	c sorties ISCCP
444
445	integer nid_isccp
446	save nid_isccp
447	cIM 090704 BEG
448	INTEGER nbapp_isccp,isccppas
449
450	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
451	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
452	DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
453
454	c cldtopres pression au sommet des nuages
455	REAL cldtopres(lmaxm1)
456	DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
457
458	INTEGER komega, nhoriRD
459
460	c taulev: numero du niveau de tau dans les sorties ISCCP
461	CHARACTER *4 taulev(kmaxm1)
462	c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
463	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
464	CHARACTER *3 pclev(lmaxm1)
465	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
466	c
467	c cnameisccp
468	CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
469	DATA cnameisccp/'pc< 50hPa, tau< 0.3',
470	. 'pc= 50-180hPa, tau< 0.3',
471	. 'pc= 180-310hPa, tau< 0.3',
472	. 'pc= 310-440hPa, tau< 0.3',
473	. 'pc= 440-560hPa, tau< 0.3',
474	. 'pc= 560-680hPa, tau< 0.3',
475	. 'pc= 680-800hPa, tau< 0.3',
476	. 'pc< 50hPa, tau= 0.3-1.3',
477	. 'pc= 50-180hPa, tau= 0.3-1.3',
478	. 'pc= 180-310hPa, tau= 0.3-1.3',
479	. 'pc= 310-440hPa, tau= 0.3-1.3',
480	. 'pc= 440-560hPa, tau= 0.3-1.3',
481	. 'pc= 560-680hPa, tau= 0.3-1.3',
482	. 'pc= 680-800hPa, tau= 0.3-1.3',
483	. 'pc< 50hPa, tau= 1.3-3.6',
484	. 'pc= 50-180hPa, tau= 1.3-3.6',
485	. 'pc= 180-310hPa, tau= 1.3-3.6',
486	. 'pc= 310-440hPa, tau= 1.3-3.6',
487	. 'pc= 440-560hPa, tau= 1.3-3.6',
488	. 'pc= 560-680hPa, tau= 1.3-3.6',
489	. 'pc= 680-800hPa, tau= 1.3-3.6',
490	. 'pc< 50hPa, tau= 3.6-9.4',
491	. 'pc= 50-180hPa, tau= 3.6-9.4',
492	. 'pc= 180-310hPa, tau= 3.6-9.4',
493	. 'pc= 310-440hPa, tau= 3.6-9.4',
494	. 'pc= 440-560hPa, tau= 3.6-9.4',
495	. 'pc= 560-680hPa, tau= 3.6-9.4',
496	. 'pc= 680-800hPa, tau= 3.6-9.4',
497	. 'pc< 50hPa, tau= 9.4-23',
498	. 'pc= 50-180hPa, tau= 9.4-23',
499	. 'pc= 180-310hPa, tau= 9.4-23',
500	. 'pc= 310-440hPa, tau= 9.4-23',
501	. 'pc= 440-560hPa, tau= 9.4-23',
502	. 'pc= 560-680hPa, tau= 9.4-23',
503	. 'pc= 680-800hPa, tau= 9.4-23',
504	. 'pc< 50hPa, tau= 23-60',
505	. 'pc= 50-180hPa, tau= 23-60',
506	. 'pc= 180-310hPa, tau= 23-60',
507	. 'pc= 310-440hPa, tau= 23-60',
508	. 'pc= 440-560hPa, tau= 23-60',
509	. 'pc= 560-680hPa, tau= 23-60',
510	. 'pc= 680-800hPa, tau= 23-60',
511	. 'pc< 50hPa, tau> 60.',
512	. 'pc= 50-180hPa, tau> 60.',
513	. 'pc= 180-310hPa, tau> 60.',
514	. 'pc= 310-440hPa, tau> 60.',
515	. 'pc= 440-560hPa, tau> 60.',
516	. 'pc= 560-680hPa, tau> 60.',
517	. 'pc= 680-800hPa, tau> 60.'/
518	c
519	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
520	c INTEGER nhorix7
521	cIM: region='3d' <==> sorties en global
522	CHARACTER*3 region
523	PARAMETER(region='3d')
524	c
525	cIM ISCCP simulator v3.4
526	c
527	logical ok_hf
528	c
529	integer nid_hf, nid_hf3d
530	save ok_hf, nid_hf, nid_hf3d
531
532	c QUESTION : noms de variables ?
533
534	#ifdef histhf
535	data ok_hf/.true./
536	#else
537	data ok_hf/.false./
538	#endif
539
540	INTEGER longcles
541	PARAMETER ( longcles = 20 )
542	REAL clesphy0( longcles )
543	c
544	c Variables quasi-arguments
545	c
546	REAL xjour
547	SAVE xjour
548	c
549	c
550	c Variables propres a la physique
551	c
552	REAL dtime
553	SAVE dtime ! pas temporel de la physique
554	c
555	INTEGER radpas
556	SAVE radpas ! frequence d'appel rayonnement
557	c
558	REAL radsol(klon)
559	SAVE radsol ! bilan radiatif au sol calcule par code radiatif
560	c
561	REAL rlat(klon)
562	SAVE rlat ! latitude pour chaque point
563	c
564	REAL rlon(klon)
565	SAVE rlon ! longitude pour chaque point
566	c
567	REAL rlonPOS(klon)
568	SAVE rlonPOS ! longitudes > 0. pour chaque point
569	c
570	cc INTEGER iflag_con
571	cc SAVE iflag_con ! indicateur de la convection
572	c
573	INTEGER itap
574	SAVE itap ! compteur pour la physique
575	c
576	REAL co2_ppm_etat0
577	c
578	REAL solaire_etat0
579	c
580	real slp(klon) ! sea level pressure
581
582	REAL ftsol(klon,nbsrf)
583	SAVE ftsol ! temperature du sol
584	c
585	REAL ftsoil(klon,nsoilmx,nbsrf)
586	SAVE ftsoil ! temperature dans le sol
587	c
588	REAL fevap(klon,nbsrf)
589	SAVE fevap ! evaporation
590	REAL fluxlat(klon,nbsrf)
591	SAVE fluxlat
592	c
593	REAL deltat(klon)
594	SAVE deltat ! ecart avec la SST de reference
595	c
596	REAL fqsurf(klon,nbsrf)
597	SAVE fqsurf ! humidite de l'air au contact de la surface
598	c
599	REAL qsol(klon)
600	SAVE qsol ! hauteur d'eau dans le sol
601	c
602	REAL fsnow(klon,nbsrf)
603	SAVE fsnow ! epaisseur neigeuse
604	c
605	REAL falbe(klon,nbsrf)
606	SAVE falbe ! albedo par type de surface
607	REAL falblw(klon,nbsrf)
608	SAVE falblw ! albedo par type de surface
609
610	c
611	c
612	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
613	c
614	REAL zmea(klon)
615	SAVE zmea ! orographie moyenne
616	c
617	REAL zstd(klon)
618	SAVE zstd ! deviation standard de l'OESM
619	c
620	REAL zsig(klon)
621	SAVE zsig ! pente de l'OESM
622	c
623	REAL zgam(klon)
624	save zgam ! anisotropie de l'OESM
625	c
626	REAL zthe(klon)
627	SAVE zthe ! orientation de l'OESM
628	c
629	REAL zpic(klon)
630	SAVE zpic ! Maximum de l'OESM
631	c
632	REAL zval(klon)
633	SAVE zval ! Minimum de l'OESM
634	c
635	REAL rugoro(klon)
636	SAVE rugoro ! longueur de rugosite de l'OESM
637	c
638	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
639	REAL zulow(klon),zvlow(klon)
640	c
641	REAL zuthe(klon),zvthe(klon)
642	SAVE zuthe
643	SAVE zvthe
644	INTEGER igwd,idx(klon),itest(klon)
645	c
646	REAL agesno(klon,nbsrf)
647	SAVE agesno ! age de la neige
648	c
649	REAL alb_neig(klon)
650	SAVE alb_neig ! albedo de la neige
651	c
652	REAL run_off_lic_0(klon)
653	SAVE run_off_lic_0
654	cKE43
655	c Variables liees a la convection de K. Emanuel (sb):
656	c
657	REAL ema_workcbmf(klon) ! cloud base mass flux
658	SAVE ema_workcbmf
659
660	REAL ema_cbmf(klon) ! cloud base mass flux
661	SAVE ema_cbmf
662
663	REAL ema_pcb(klon) ! cloud base pressure
664	SAVE ema_pcb
665
666	REAL ema_pct(klon) ! cloud top pressure
667	SAVE ema_pct
668
669	REAL bas, top ! cloud base and top levels
670	SAVE bas
671	SAVE top
672
673	REAL Ma(klon,klev) ! undilute upward mass flux
674	SAVE Ma
675	REAL qcondc(klon,klev) ! in-cld water content from convect
676	SAVE qcondc
677	REAL ema_work1(klon, klev), ema_work2(klon, klev)
678	SAVE ema_work1, ema_work2
679	REAL wdn(klon), tdn(klon), qdn(klon)
680
681	REAL wd(klon) ! sb
682	SAVE wd ! sb
683
684	c Variables locales pour la couche limite (al1):
685	c
686	cAl1 REAL pblh(klon) ! Hauteur de couche limite
687	cAl1 SAVE pblh
688	c34EK
689	c
690	c Variables locales:
691	c
692	REAL cdragh(klon) ! drag coefficient pour T and Q
693	REAL cdragm(klon) ! drag coefficient pour vent
694	cAA
695	cAA Pour phytrac
696	cAA
697	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
698	REAL yu1(klon) ! vents dans la premiere couche U
699	REAL yv1(klon) ! vents dans la premiere couche V
700	REAL ffonte(klon,nbsrf) !Flux thermique utilise pour fondre la neige
701	REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface
702	c !et necessaire pour limiter la
703	c !hauteur de neige, en kg/m2/s
704	REAL zxffonte(klon), zxfqcalving(klon)
705
706	c$$$ LOGICAL offline ! Controle du stockage ds "physique"
707	c$$$ PARAMETER (offline=.false.)
708	c$$$ INTEGER physid
709	REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
710	save pfrac_impa
711	REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
712	save pfrac_nucl
713	REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
714	save pfrac_1nucl
715	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
716	REAL frac_nucl(klon,klev) ! idem (nucleation)
717	#ifdef INCA
718	INTEGER :: iii
719	REAL :: calday
720	#endif
721
722	cAA
723	REAL rain_fall(klon) ! pluie
724	REAL snow_fall(klon) ! neige
725	save snow_fall, rain_fall
726	cIM cf FH pour Tiedtke 080604
727	REAL rain_tiedtke(klon),snow_tiedtke(klon)
728	c
729	REAL total_rain(klon), nday_rain(klon)
730	save nday_rain
731	c
732	REAL evap(klon), devap(klon) ! evaporation et sa derivee
733	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
734	REAL dlw(klon) ! derivee infra rouge
735	cym
736	SAVE dlw
737	cym
738	REAL bils(klon) ! bilan de chaleur au sol
739	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
740	C ! type de sous-surface et pondere par la fraction
741	REAL fder(klon) ! Derive de flux (sensible et latente)
742	save fder
743	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
744	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
745	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
746	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
747	c
748	REAL frugs(klon,nbsrf) ! longueur de rugosite
749	save frugs
750	REAL zxrugs(klon) ! longueur de rugosite
751	c
752	c Conditions aux limites
753	c
754	INTEGER julien
755	c
756	INTEGER lmt_pas
757	SAVE lmt_pas ! frequence de mise a jour
758	REAL pctsrf(klon,nbsrf)
759	cIM
760	REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
761	REAL paire_ter(klon) !surfaces terre
762	c
763	SAVE pctsrf ! sous-fraction du sol
764	REAL albsol(klon)
765	SAVE albsol ! albedo du sol total
766	REAL albsollw(klon)
767	SAVE albsollw ! albedo du sol total
768
769	REAL wo(klon,klev)
770	SAVE wo ! ozone
771	c======================================================================
772	c
773	c Declaration des procedures appelees
774	c
775	EXTERNAL angle ! calculer angle zenithal du soleil
776	EXTERNAL alboc ! calculer l'albedo sur ocean
777	EXTERNAL ajsec ! ajustement sec
778	EXTERNAL clmain ! couche limite
779	EXTERNAL conlmd ! convection (schema LMD)
780	cKE43
781	EXTERNAL conema3 ! convect4.3
782	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
783	cAA
784	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
785	c ! stockage des coefficients necessaires au
786	c ! lessivage OFF-LINE et ON-LINE
787	EXTERNAL hgardfou ! verifier les temperatures
788	EXTERNAL nuage ! calculer les proprietes radiatives
789	EXTERNAL o3cm ! initialiser l'ozone
790	EXTERNAL orbite ! calculer l'orbite terrestre
791	EXTERNAL ozonecm ! prescrire l'ozone
792	EXTERNAL phyetat0 ! lire l'etat initial de la physique
793	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
794	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
795	EXTERNAL suphec ! initialiser certaines constantes
796	EXTERNAL transp ! transport total de l'eau et de l'energie
797	EXTERNAL ecribina ! ecrire le fichier binaire global
798	EXTERNAL ecribins ! ecrire le fichier binaire global
799	EXTERNAL ecrirega ! ecrire le fichier binaire regional
800	EXTERNAL ecriregs ! ecrire le fichier binaire regional
801	cIM
802	EXTERNAL haut2bas !variables de haut en bas
803	INTEGER lnblnk1
804	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
805	!caracter
806	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
807	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
808	c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
809	c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
810	c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
811	c
812	c Variables locales
813	c
814	real clwcon(klon,klev),rnebcon(klon,klev)
815	real clwcon0(klon,klev),rnebcon0(klon,klev)
816	cIM cf. AM 081204 BEG
817	real clwcon0th(klon,klev),rnebcon0th(klon,klev)
818	cIM cf. AM 081204 END
819	save rnebcon, clwcon
820
821	REAL rhcl(klon,klev) ! humiditi relative ciel clair
822	REAL dialiq(klon,klev) ! eau liquide nuageuse
823	REAL diafra(klon,klev) ! fraction nuageuse
824	REAL cldliq(klon,klev) ! eau liquide nuageuse
825	REAL cldfra(klon,klev) ! fraction nuageuse
826	REAL cldtau(klon,klev) ! epaisseur optique
827	REAL cldemi(klon,klev) ! emissivite infrarouge
828	c
829	CXXX PB
830	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
831	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
832	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
833	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
834	c
835	REAL zxfluxt(klon, klev)
836	REAL zxfluxq(klon, klev)
837	REAL zxfluxu(klon, klev)
838	REAL zxfluxv(klon, klev)
839	CXXX
840	REAL heat(klon,klev) ! chauffage solaire
841	REAL heat0(klon,klev) ! chauffage solaire ciel clair
842	REAL cool(klon,klev) ! refroidissement infrarouge
843	REAL cool0(klon,klev) ! refroidissement infrarouge ciel clair
844	REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
845	real sollwdown(klon) ! downward LW flux at surface
846	cIM BEG
847	real sollwdownclr(klon) ! downward CS LW flux at surface
848	real toplwdown(klon) ! downward CS LW flux at TOA
849	real toplwdownclr(klon) ! downward CS LW flux at TOA
850	cIM END
851	REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
852	REAL albpla(klon)
853	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
854	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
855	c Le rayonnement n'est pas calcule tous les pas, il faut donc
856	c sauvegarder les sorties du rayonnement
857	SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
858	SAVE sollwdownclr, toplwdown, toplwdownclr
859	SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
860	c
861	INTEGER itaprad
862	SAVE itaprad
863	c
864	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
865	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
866	c
867	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
868	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
869	c
870	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
871	c
872	REAL dist, rmu0(klon), fract(klon)
873	REAL zdtime, zlongi
874	c
875	CHARACTER*2 str2
876	CHARACTER*2 iqn
877	c
878	REAL qcheck
879	REAL z_avant(klon), z_apres(klon), z_factor(klon)
880	LOGICAL zx_ajustq
881	c
882	REAL za, zb
883	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
884	real zqsat(klon,klev)
885	INTEGER i, k, iq, ig, j, nsrf, ll, iiq
886	REAL t_coup
887	PARAMETER (t_coup=234.0)
888	c
889	REAL zphi(klon,klev)
890	REAL zx_relief(iim,jjmp1)
891	REAL zx_aire(iim,jjmp1)
892	c
893	cIM cf. AM Variables locales pour la CLA (hbtm2)
894	c
895	REAL pblh(klon, nbsrf) ! Hauteur de couche limite
896	REAL plcl(klon, nbsrf) ! Niveau de condensation de la CLA
897	REAL capCL(klon, nbsrf) ! CAPE de couche limite
898	REAL oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
899	REAL cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
900	REAL pblt(klon, nbsrf) ! T a la Hauteur de couche limite
901	REAL therm(klon, nbsrf)
902	REAL trmb1(klon, nbsrf) ! deep_cape
903	REAL trmb2(klon, nbsrf) ! inhibition
904	REAL trmb3(klon, nbsrf) ! Point Omega
905	c Grdeurs de sorties
906	REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
907	REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
908	REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
909	REAL s_trmb3(klon)
910	cKE43
911	c Variables locales pour la convection de K. Emanuel (sb):
912	c
913	REAL upwd(klon,klev) ! saturated updraft mass flux
914	REAL dnwd(klon,klev) ! saturated downdraft mass flux
915	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
916	REAL tvp(klon,klev) ! virtual temp of lifted parcel
917	REAL cape(klon) ! CAPE
918	SAVE cape
919	CHARACTER*40 capemaxcels !max(CAPE)
920
921	REAL pbase(klon) ! cloud base pressure
922	SAVE pbase
923	REAL bbase(klon) ! cloud base buoyancy
924	SAVE bbase
925	REAL rflag(klon) ! flag fonctionnement de convect
926	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
927	c -- convect43:
928	INTEGER ntra ! nb traceurs pour convect4.3
929	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
930	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
931	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
932	REAL dplcldt(klon), dplcldr(klon)
933	c? . condm_con(klon,klev),conda_con(klon,klev),
934	c? . mr_con(klon,klev),ep_con(klon,klev)
935	c? . ,sadiab(klon,klev),wadiab(klon,klev)
936	c --
937	c34EK
938	c
939	c Variables du changement
940	c
941	c con: convection
942	c lsc: condensation a grande echelle (Large-Scale-Condensation)
943	c ajs: ajustement sec
944	c eva: evaporation de l'eau liquide nuageuse
945	c vdf: couche limite (Vertical DiFfusion)
946	REAL d_t_con(klon,klev),d_q_con(klon,klev)
947	REAL d_u_con(klon,klev),d_v_con(klon,klev)
948	REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
949	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
950	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
951	REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
952	REAL d_t_oli(klon,klev) !tendances dues a oro et lif
953	REAL rneb(klon,klev)
954	c
955	*********************************************************
956	* declarations
957	real zqasc(klon,klev)
958	save zqasc
959
960	*********************************************************
961	cIM 081204 END
962	c
963	REAL pmfu(klon,klev), pmfd(klon,klev)
964	REAL pen_u(klon,klev), pen_d(klon,klev)
965	REAL pde_u(klon,klev), pde_d(klon,klev)
966	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
967	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
968	REAL prfl(klon,klev+1), psfl(klon,klev+1)
969	c
970	INTEGER ibas_con(klon), itop_con(klon)
971	cym
972	SAVE ibas_con,itop_con
973	cym
974	REAL rain_con(klon), rain_lsc(klon)
975	REAL snow_con(klon), snow_lsc(klon)
976	REAL d_ts(klon,nbsrf)
977	c
978	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
979	REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
980	c
981	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
982	REAL d_t_oro(klon,klev)
983	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
984	REAL d_t_lif(klon,klev)
985	REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif
986
987	REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
988	real ratqsbas,ratqshaut
989	save ratqsbas,ratqshaut, ratqs
990	real zpt_conv(klon,klev)
991
992	c Parametres lies au nouveau schema de nuages (SB, PDF)
993	real fact_cldcon
994	real facttemps
995	logical ok_newmicro
996	save ok_newmicro
997	save fact_cldcon,facttemps
998	real facteur
999
1000	integer iflag_cldcon
1001	save iflag_cldcon
1002
1003	logical ptconv(klon,klev)
1004	cIM cf. AM 081204 BEG
1005	logical ptconvth(klon,klev)
1006	cIM cf. AM 081204 END
1007	c
1008	c Variables liees a l'ecriture de la bande histoire physique
1009	c
1010	c======================================================================
1011	c
1012	cIM cf. AM 081204 BEG
1013	c declarations pour sortir sur une sous-region
1014	integer imin_ins,imax_ins,jmin_ins,jmax_ins
1015	save imin_ins,imax_ins,jmin_ins,jmax_ins
1016	c real lonmin_ins,lonmax_ins,latmin_ins
1017	c s ,latmax_ins
1018	c data lonmin_ins,lonmax_ins,latmin_ins
1019	c s ,latmax_ins/
1020	c valeurs initiales s -5.,20.,41.,55./
1021	c s 100.,130.,-20.,20./
1022	c s -180.,180.,-90.,90./
1023	c======================================================================
1024	cIM cf. AM 081204 END
1025
1026	c
1027	integer itau_w ! pas de temps ecriture = itap + itau_phy
1028	c
1029	c
1030	c Variables locales pour effectuer les appels en serie
1031	c
1032	REAL t_seri(klon,klev), q_seri(klon,klev)
1033	REAL ql_seri(klon,klev),qs_seri(klon,klev)
1034	REAL u_seri(klon,klev), v_seri(klon,klev)
1035	c
1036	REAL tr_seri(klon,klev,nbtr)
1037	REAL d_tr(klon,klev,nbtr)
1038
1039	REAL zx_rh(klon,klev)
1040
1041	INTEGER length
1042	PARAMETER ( length = 100 )
1043	REAL tabcntr0( length )
1044	c
1045	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
1046	c
1047	cIM AMIP2 BEG
1048	REAL moyglo, mountor
1049	cIM 141004 BEG
1050	REAL zustrdr(klon), zvstrdr(klon)
1051	REAL zustrli(klon), zvstrli(klon)
1052	REAL zustrph(klon), zvstrph(klon)
1053	REAL aam, torsfc
1054	cIM 141004 END
1055	cIM 190504 BEG
1056	INTEGER ij, imp1jmp1
1057	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
1058	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
1059	REAL padyn(iim+1,jjmp1,klev+1)
1060	REAL dudyn(iim+1,jjmp1,klev)
1061	REAL rlatdyn(iim+1,jjmp1)
1062	cIM 190504 END
1063	LOGICAL ok_msk
1064	REAL msk(klon)
1065	cIM
1066	REAL airetot, pi
1067	REAL zm_wo(jjmp1, klev)
1068	cIM AMIP2 END
1069	c
1070	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
1071	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
1072	REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
1073	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
1074	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
1075	c
1076	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
1077	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
1078	c
1079	cIM 280405 BEG
1080	INTEGER nid_bilKPins, nid_bilKPave
1081	SAVE nid_bilKPins, nid_bilKPave
1082	c
1083	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1084	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1085	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1086	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1087	c
1088	cIM 280405 END
1089	c
1090	INTEGER nhori, nvert, nvert1
1091	c REAL zstok
1092	REAL zsto, zout, zsto1, zsto2
1093	c REAL zstoave, zstoin
1094	REAL zstophy, zstorad, zstohf, zstoday, zstomth
1095	real zjulian
1096	save zjulian
1097
1098	character*20 modname
1099	character*80 abort_message
1100	logical ok_sync
1101	real date0
1102	integer idayref
1103
1104	C essai writephys
1105	integer fid_day, fid_mth, fid_ins
1106	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1107	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1108	parameter (prof2d_on = 1, prof3d_on = 2,
1109	. prof2d_av = 3, prof3d_av = 4)
1110	character*30 nom_fichier
1111	character*10 varname
1112	character*40 vartitle
1113	character*20 varunits
1114	C Variables liees au bilan d'energie et d'enthalpi
1115	REAL ztsol(klon)
1116	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1117	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1118	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1119	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1120	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1121	REAL d_h_vcol_phy
1122	REAL fs_bound, fq_bound
1123	SAVE d_h_vcol_phy
1124	REAL zero_v(klon)
1125	CHARACTER*15 ztit
1126	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1127	SAVE ip_ebil
1128	DATA ip_ebil/0/
1129	INTEGER if_ebil ! level for energy conserv. dignostics
1130	SAVE if_ebil
1131	c+jld ec_conser
1132	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
1133	REAL ZRCPD
1134	c-jld ec_conser
1135	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1136	REAL t2m(klon,nbsrf), q2m(klon,nbsrf) !temperature, humidite a 2m
1137	REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m
1138	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1139	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1140	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1141	CHARACTER*40 tinst, tave, typeval
1142	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
1143	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
1144	REAL sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
1145	SAVE sulfate_pi
1146
1147	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1148
1149	REAL re(klon, klev) ! Cloud droplet effective radius
1150	REAL fl(klon, klev) ! denominator of re
1151
1152	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1153
1154	! Aerosol optical properties
1155	REAL tau_ae(klon,klev,2), piz_ae(klon,klev,2)
1156	REAL cg_ae(klon,klev,2)
1157
1158	REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
1159	! ok_ade=T -ADE=topswad-topsw
1160
1161	REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
1162	! ok_aie=T ->
1163	! ok_ade=T -AIE=topswai-topswad
1164	! ok_ade=F -AIE=topswai-topsw
1165
1166	REAL aerindex(klon) ! POLDER aerosol index
1167
1168	! Parameters
1169	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1170	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1171	cym
1172	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1173	cym
1174	cjq-end
1175	c
1176	c Declaration des constantes et des fonctions thermodynamiques
1177	c
1178	#include "YOMCST.h"
1179	#include "YOETHF.h"
1180	#include "FCTTRE.h"
1181	c======================================================================
1182	modname = 'physiq'
1183	IF (if_ebil.ge.1) THEN
1184	DO i=1,klon
1185	zero_v(i)=0.
1186	END DO
1187	END IF
1188	ok_sync=.TRUE.
1189	IF (nqmax .LT. 2) THEN
1190	abort_message = 'eaux vapeur et liquide sont indispensables'
1191	CALL abort_gcm (modname,abort_message,1)
1192	ENDIF
1193	IF (debut) THEN
1194	CALL suphec ! initialiser constantes et parametres phys.
1195	ENDIF
1196
1197
1198	c======================================================================
1199	xjour = rjourvrai
1200	c
1201	c Si c'est le debut, il faut initialiser plusieurs choses
1202	c ********
1203	c
1204	IF (debut) THEN
1205	C
1206	!rv
1207	u10m(:,:)=0.
1208	v10m(:,:)=0.
1209	t2m(:,:)=0.
1210	q2m(:,:)=0.
1211	ffonte(:,:)=0.
1212	fqcalving(:,:)=0.
1213	piz_ae(:,:,:)=0.
1214	tau_ae(:,:,:)=0.
1215	cg_ae(:,:,:)=0.
1216	rain_con(:)=0.
1217	snow_con(:)=0.
1218	bl95_b0=0.
1219	bl95_b1=0.
1220	topswai(:)=0.
1221	topswad(:)=0.
1222	solswai(:)=0.
1223	solswad(:)=0.
1224	crv
1225	c anne
1226	d_u_con(:,:) = 0.0
1227	d_v_con(:,:) = 0.0
1228	rnebcon0(:,:) = 0.0
1229	clwcon0(:,:) = 0.0
1230	rnebcon(:,:) = 0.0
1231	clwcon(:,:) = 0.0
1232	paire_ter(:) = 0.0
1233	c fin anne
1234
1235	cym
1236	wfbils(:,:)=0
1237	cym
1238	IF (if_ebil.ge.1) d_h_vcol_phy=0.
1239	c
1240	c appel a la lecture du run.def physique
1241	c
1242	call conf_phys(pdtphys, ocean, ok_veget, ok_journe, ok_mensuel,
1243	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
1244	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
1245	. ok_ade, ok_aie,
1246	. bl95_b0, bl95_b1,
1247	. iflag_thermals,nsplit_thermals)
1248
1249	c
1250	c
1251	c Initialiser les compteurs:
1252	c
1253
1254	frugs = 0.
1255	itap = 0
1256	itaprad = 0
1257	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1258	. rlat,rlon,pctsrf, ftsol,ftsoil,
1259	cIM "slab" ocean
1260	. ocean, tslab,seaice,
1261	. fqsurf,qsol,fsnow,
1262	. falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
1263	. dlw,radsol,frugs,agesno,clesphy0,
1264	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1265	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1266	. run_off_lic_0)
1267
1268	c ATTENTION : il faudra a terme relire q2 dans l'etat initial
1269	q2(:,:,:)=1.e-8
1270	c
1271	radpas = NINT( 86400./dtime/nbapp_rad)
1272	c
1273	C on remet le calendrier a zero
1274	c
1275	IF (raz_date .eq. 1) THEN
1276	itau_phy = 0
1277	ENDIF
1278
1279	cIM cf. AM 081204 BEG
1280	PRINT*,'cycle_diurne3 =',cycle_diurne
1281	cIM cf. AM 081204 END
1282	c
1283	IF(ocean.NE.'force ') THEN
1284	ok_ocean=.TRUE.
1285	ENDIF
1286	c
1287	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
1288	, ok_instan, ok_region )
1289	c
1290	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1291	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1292	. pdtphys
1293	abort_message='Pas physique n est pas correct '
1294	call abort_gcm(modname,abort_message,1)
1295	ENDIF
1296	IF (nlon .NE. klon) THEN
1297	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1298	. klon
1299	abort_message='nlon et klon ne sont pas coherents'
1300	call abort_gcm(modname,abort_message,1)
1301	ENDIF
1302	IF (nlev .NE. klev) THEN
1303	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1304	. klev
1305	abort_message='nlev et klev ne sont pas coherents'
1306	call abort_gcm(modname,abort_message,1)
1307	ENDIF
1308	c
1309	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1310	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1311	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1312	abort_message='Nbre d appels au rayonnement insuffisant'
1313	call abort_gcm(modname,abort_message,1)
1314	ENDIF
1315	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1316	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1317	. ok_cvl
1318	c
1319	cKE43
1320	c Initialisation pour la convection de K.E. (sb):
1321	IF (iflag_con.GE.3) THEN
1322
1323	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1324
1325	DO i = 1, klon
1326	ema_cbmf(i) = 0.
1327	ema_pcb(i) = 0.
1328	ema_pct(i) = 0.
1329	ema_workcbmf(i) = 0.
1330	ENDDO
1331	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1332	DO i = 1, klon
1333	ibas_con(i) = 1
1334	itop_con(i) = 1
1335	ENDDO
1336	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1337
1338	ENDIF
1339
1340	c34EK
1341	IF (ok_orodr) THEN
1342	DO i=1,klon
1343	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1344	ENDDO
1345	CALL SUGWD(klon,klev,paprs,pplay)
1346	DO i=1,klon
1347	zuthe(i)=0.
1348	zvthe(i)=0.
1349	if(zstd(i).gt.10.)then
1350	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1351	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1352	endif
1353	ENDDO
1354	ENDIF
1355	c
1356	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1357	WRITE(lunout,*)'La frequence de lecture surface est de ',
1358	. lmt_pas
1359	c
1360	c Initialiser le couplage si necessaire
1361	c
1362	npas = 0
1363	nexca = 0
1364	if (ocean == 'couple') then
1365	npas = itaufin/ iphysiq
1366	nexca = 86400 / dtime
1367	write(lunout,*)' ##### Ocean couple #####'
1368	write(lunout,*)' Valeurs des pas de temps'
1369	write(lunout,*)' npas = ', npas
1370	write(lunout,*)' nexca = ', nexca
1371	endif
1372	c
1373	capemaxcels = 't_max(X)'
1374	t2mincels = 't_min(X)'
1375	t2maxcels = 't_max(X)'
1376	tinst = 'inst(X)'
1377	tave = 'ave(X)'
1378	cIM cf. AM 081204 BEG
1379	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1380	cIM cf. AM 081204 END
1381	c
1382	c=============================================================
1383	c Initialisation des sorties
1384	c=============================================================
1385
1386	#ifdef CPP_IOIPSL
1387
1388	#ifdef histhf
1389	#include "ini_histhf.h"
1390	#endif
1391
1392	#ifdef histday
1393	#include "ini_histday.h"
1394	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1395	c#include "ini_bilKP_ins.h"
1396	c#include "ini_bilKP_ave.h"
1397	#include "ini_histday_seri.h"
1398	#endif
1399
1400	#ifdef histmth
1401	#include "ini_histmth.h"
1402	#endif
1403
1404	#ifdef histins
1405	#include "ini_histins.h"
1406	#endif
1407
1408	#ifdef histISCCP
1409	#include "ini_histISCCP.h"
1410	#endif
1411
1412	#ifdef histmthNMC
1413	#include "ini_histmthNMC.h"
1414	#endif
1415
1416	#endif
1417
1418	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1419	date0 = zjulian
1420	C date0 = day_ini
1421	WRITE(,) 'physiq date0 : ',date0
1422	c
1423	c
1424	c
1425	c Prescrire l'ozone dans l'atmosphere
1426	c
1427	c
1428	cc DO i = 1, klon
1429	cc DO k = 1, klev
1430	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1431	cc ENDDO
1432	cc ENDDO
1433	c
1434	#ifdef INCA
1435	iii = MOD(NINT(xjour),360)
1436	calday = FLOAT(iii) + gmtime
1437	WRITE(lunout,*) 'initial time ', xjour, calday
1438	#ifdef INCAINFO
1439	WRITE(lunout,*) 'Appel CHEMINI ...'
1440	#endif
1441	CALL chemini( rpi,
1442	$ rg,
1443	$ ra,
1444	$ airephy,
1445	$ rlat,
1446	$ rlon,
1447	$ presnivs,
1448	$ calday,
1449	$ klon,
1450	$ nqmax,
1451	$ pdtphys,
1452	$ annee_ref,
1453	$ day_ini)
1454	#ifdef INCAINFO
1455	WRITE(lunout,*) 'OK.'
1456	#endif
1457	#endif
1458	c
1459	ENDIF
1460	c
1461	c ************** Fin de IF ( debut ) *************
1462	c
1463	c
1464	c Mettre a zero des variables de sortie (pour securite)
1465	c
1466	DO i = 1, klon
1467	d_ps(i) = 0.0
1468	ENDDO
1469	DO k = 1, klev
1470	DO i = 1, klon
1471	d_t(i,k) = 0.0
1472	d_u(i,k) = 0.0
1473	d_v(i,k) = 0.0
1474	ENDDO
1475	ENDDO
1476	DO iq = 1, nqmax
1477	DO k = 1, klev
1478	DO i = 1, klon
1479	d_qx(i,k,iq) = 0.0
1480	ENDDO
1481	ENDDO
1482	ENDDO
1483	c
1484	c Ne pas affecter les valeurs entrees de u, v, h, et q
1485	c
1486	DO k = 1, klev
1487	DO i = 1, klon
1488	t_seri(i,k) = t(i,k)
1489	u_seri(i,k) = u(i,k)
1490	v_seri(i,k) = v(i,k)
1491	q_seri(i,k) = qx(i,k,ivap)
1492	ql_seri(i,k) = qx(i,k,iliq)
1493	qs_seri(i,k) = 0.
1494	ENDDO
1495	ENDDO
1496	IF (nqmax.GE.3) THEN
1497	DO iq = 3, nqmax
1498	DO k = 1, klev
1499	DO i = 1, klon
1500	tr_seri(i,k,iq-2) = qx(i,k,iq)
1501	ENDDO
1502	ENDDO
1503	ENDDO
1504	ELSE
1505	DO k = 1, klev
1506	DO i = 1, klon
1507	tr_seri(i,k,1) = 0.0
1508	ENDDO
1509	ENDDO
1510	ENDIF
1511	C
1512	DO i = 1, klon
1513	ztsol(i) = 0.
1514	ENDDO
1515	DO nsrf = 1, nbsrf
1516	DO i = 1, klon
1517	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1518	ENDDO
1519	ENDDO
1520	C
1521	IF (if_ebil.ge.1) THEN
1522	ztit='after dynamic'
1523	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1524	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1525	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1526	C Comme les tendances de la physique sont ajoute dans la dynamique,
1527	C on devrait avoir que la variation d'entalpie par la dynamique
1528	C est egale a la variation de la physique au pas de temps precedent.
1529	C Donc la somme de ces 2 variations devrait etre nulle.
1530	call diagphy(airephy,ztit,ip_ebil
1531	e , zero_v, zero_v, zero_v, zero_v, zero_v
1532	e , zero_v, zero_v, zero_v, ztsol
1533	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1534	s , fs_bound, fq_bound )
1535	END IF
1536
1537	c Diagnostiquer la tendance dynamique
1538	c
1539	IF (ancien_ok) THEN
1540	DO k = 1, klev
1541	DO i = 1, klon
1542	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1543	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1544	ENDDO
1545	ENDDO
1546	ELSE
1547	DO k = 1, klev
1548	DO i = 1, klon
1549	d_t_dyn(i,k) = 0.0
1550	d_q_dyn(i,k) = 0.0
1551	ENDDO
1552	ENDDO
1553	ancien_ok = .TRUE.
1554	ENDIF
1555	c
1556	c Ajouter le geopotentiel du sol:
1557	c
1558	DO k = 1, klev
1559	DO i = 1, klon
1560	zphi(i,k) = pphi(i,k) + pphis(i)
1561	ENDDO
1562	ENDDO
1563	c
1564	c Verifier les temperatures
1565	c
1566	CALL hgardfou(t_seri,ftsol,'debutphy')
1567	c
1568	c Incrementer le compteur de la physique
1569	c
1570	itap = itap + 1
1571	julien = MOD(NINT(xjour),360)
1572	if (julien .eq. 0) julien = 360
1573	c
1574	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1575	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1576	c
1577	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1578	WRITE(lunout,*)' PHYS cond julien ',julien
1579	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1580	ENDIF
1581	c
1582	c Re-evaporer l'eau liquide nuageuse
1583	c
1584	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1585	DO i = 1, klon
1586	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1587	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1588	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1589	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1590	zb = MAX(0.0,ql_seri(i,k))
1591	za = - MAX(0.0,ql_seri(i,k))
1592	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1593	t_seri(i,k) = t_seri(i,k) + za
1594	q_seri(i,k) = q_seri(i,k) + zb
1595	ql_seri(i,k) = 0.0
1596	d_t_eva(i,k) = za
1597	d_q_eva(i,k) = zb
1598	ENDDO
1599	ENDDO
1600	c
1601	IF (if_ebil.ge.2) THEN
1602	ztit='after reevap'
1603	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
1604	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1605	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1606	call diagphy(airephy,ztit,ip_ebil
1607	e , zero_v, zero_v, zero_v, zero_v, zero_v
1608	e , zero_v, zero_v, zero_v, ztsol
1609	e , d_h_vcol, d_qt, d_ec
1610	s , fs_bound, fq_bound )
1611	C
1612	END IF
1613	C
1614	c
1615	c Appeler la diffusion verticale (programme de couche limite)
1616	c
1617	DO i = 1, klon
1618	c if (.not. ok_veget) then
1619	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1620	c endif
1621	c frugs(i,is_lic) = rugoro(i)
1622	c frugs(i,is_oce) = rugmer(i)
1623	c frugs(i,is_sic) = 0.001
1624	zxrugs(i) = 0.0
1625	ENDDO
1626	DO nsrf = 1, nbsrf
1627	DO i = 1, klon
1628	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1629	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1630	ENDDO
1631	ENDDO
1632	DO nsrf = 1, nbsrf
1633	DO i = 1, klon
1634	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1635	ENDDO
1636	ENDDO
1637	c
1638	C calculs necessaires au calcul de l'albedo dans l'interface
1639	c
1640	CALL orbite(FLOAT(julien),zlongi,dist)
1641	IF (cycle_diurne) THEN
1642	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1643	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1644	ELSE
1645	rmu0 = -999.999
1646	ENDIF
1647	c
1648	C Calcul de l'abedo moyen par maille
1649	albsol(:)=0.
1650	albsollw(:)=0.
1651	DO nsrf = 1, nbsrf
1652	DO i = 1, klon
1653	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1654	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1655	ENDDO
1656	ENDDO
1657	C
1658	C Repartition sous maille des flux LW et SW
1659	C Modif OM+PASB+JLD
1660	C Repartition du longwave par sous-surface linearisee
1661	Cn
1662
1663	DO nsrf = 1, nbsrf
1664	DO i = 1, klon
1665	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1666	c$$$ fsollw(i,nsrf) = sollw(i)
1667	fsollw(i,nsrf) = sollw(i)
1668	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
1669	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1670	ENDDO
1671	ENDDO
1672
1673	fder = dlw
1674
1675
1676	CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
1677	e t_seri,q_seri,u_seri,v_seri,
1678	e julien, rmu0, co2_ppm,
1679	e ok_veget, ocean, npas, nexca, ftsol,
1680	$ soil_model,cdmmax, cdhmax,
1681	$ ksta, ksta_ter, ok_kzmin, ftsoil, qsol,
1682	cIM BAD $ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1683	$ paprs,pplay, fsnow,fqsurf,fevap,falbe,falblw,
1684	$ fluxlat,
1685	e rain_fall, snow_fall,
1686	e fsolsw, fsollw, sollwdown, fder,
1687	e rlon, rlat, cuphy, cvphy, frugs,
1688	e debut, lafin, agesno,rugoro ,
1689	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1690	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1691	s q2,
1692	s dsens, devap,
1693	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1694	cIM cf. AM 081204 BEG
1695	s pblh,capCL,oliqCL,cteiCL,pblT,
1696	s therm,trmb1,trmb2,trmb3,plcl,
1697	cIM cf. AM 081204 END
1698	s fqcalving, ffonte, run_off_lic_0,
1699	cIM "slab" ocean
1700	s fluxo, fluxg, tslab, seaice)
1701	c
1702	CXXX PB
1703	CXXX Incrementation des flux
1704	CXXX
1705
1706	zxfluxt=0.
1707	zxfluxq=0.
1708	zxfluxu=0.
1709	zxfluxv=0.
1710	DO nsrf = 1, nbsrf
1711	DO k = 1, klev
1712	DO i = 1, klon
1713	zxfluxt(i,k) = zxfluxt(i,k) +
1714	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1715	zxfluxq(i,k) = zxfluxq(i,k) +
1716	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1717	zxfluxu(i,k) = zxfluxu(i,k) +
1718	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1719	zxfluxv(i,k) = zxfluxv(i,k) +
1720	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1721	END DO
1722	END DO
1723	END DO
1724	DO i = 1, klon
1725	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1726	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1727	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1728	fder(i) = dlw(i) + dsens(i) + devap(i)
1729	ENDDO
1730
1731
1732	DO k = 1, klev
1733	DO i = 1, klon
1734	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1735	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1736	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1737	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1738	ENDDO
1739	ENDDO
1740	c
1741	IF (if_ebil.ge.2) THEN
1742	ztit='after clmain'
1743	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1744	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1745	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1746	call diagphy(airephy,ztit,ip_ebil
1747	e , zero_v, zero_v, zero_v, zero_v, sens
1748	e , evap , zero_v, zero_v, ztsol
1749	e , d_h_vcol, d_qt, d_ec
1750	s , fs_bound, fq_bound )
1751	END IF
1752	C
1753	c
1754	c Incrementer la temperature du sol
1755	c
1756	DO i = 1, klon
1757	zxtsol(i) = 0.0
1758	zxfluxlat(i) = 0.0
1759	c
1760	zt2m(i) = 0.0
1761	zq2m(i) = 0.0
1762	zu10m(i) = 0.0
1763	zv10m(i) = 0.0
1764	cIM cf JLD ??
1765	zxffonte(i) = 0.0
1766	zxfqcalving(i) = 0.0
1767	cIM cf. AM 081204 BEG
1768	c
1769	s_pblh(i) = 0.0
1770	s_lcl(i) = 0.0
1771	s_capCL(i) = 0.0
1772	s_oliqCL(i) = 0.0
1773	s_cteiCL(i) = 0.0
1774	s_pblT(i) = 0.0
1775	s_therm(i) = 0.0
1776	s_trmb1(i) = 0.0
1777	s_trmb2(i) = 0.0
1778	s_trmb3(i) = 0.0
1779	c
1780	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1781	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1782	$ THEN
1783	WRITE(,) 'physiq : pb sous surface au point ', i,
1784	$ pctsrf(i, 1 : nbsrf)
1785	ENDIF
1786	ENDDO
1787	DO nsrf = 1, nbsrf
1788	DO i = 1, klon
1789	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1790	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1791	cIM cf. JLD
1792	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1793	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1794	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1795	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1796	cccIM
1797	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1798	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1799	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1800	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1801	cIM cf JLD ??
1802	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
1803	zxfqcalving(i) = zxfqcalving(i) +
1804	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
1805	cIM cf. AM 081204 BEG
1806	s_pblh(i) = s_pblh(i) + pblh(i,nsrf)*pctsrf(i,nsrf)
1807	s_lcl(i) = s_lcl(i) + plcl(i,nsrf)*pctsrf(i,nsrf)
1808	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) *pctsrf(i,nsrf)
1809	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf) *pctsrf(i,nsrf)
1810	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf) *pctsrf(i,nsrf)
1811	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) *pctsrf(i,nsrf)
1812	s_therm(i) = s_therm(i) + therm(i,nsrf) *pctsrf(i,nsrf)
1813	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) *pctsrf(i,nsrf)
1814	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) *pctsrf(i,nsrf)
1815	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) *pctsrf(i,nsrf)
1816	c ENDIF
1817	ENDDO
1818	ENDDO
1819
1820	c
1821	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1822	c
1823	DO nsrf = 1, nbsrf
1824	DO i = 1, klon
1825	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1826	cccIM
1827	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1828	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1829	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1830	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1831	cIM cf JLD ??
1832	IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
1833	IF (pctsrf(i,nsrf) .LT. epsfra)
1834	. fqcalving(i,nsrf) = zxfqcalving(i)
1835	cIM cf. AM 081204 BEG
1836	IF (pctsrf(i,nsrf) .LT. epsfra) pblh(i,nsrf)=s_pblh(i)
1837	IF (pctsrf(i,nsrf) .LT. epsfra) plcl(i,nsrf)=s_lcl(i)
1838	IF (pctsrf(i,nsrf) .LT. epsfra) capCL(i,nsrf)=s_capCL(i)
1839	IF (pctsrf(i,nsrf) .LT. epsfra) oliqCL(i,nsrf)=s_oliqCL(i)
1840	IF (pctsrf(i,nsrf) .LT. epsfra) cteiCL(i,nsrf)=s_cteiCL(i)
1841	IF (pctsrf(i,nsrf) .LT. epsfra) pblT(i,nsrf)=s_pblT(i)
1842	IF (pctsrf(i,nsrf) .LT. epsfra) therm(i,nsrf)=s_therm(i)
1843	IF (pctsrf(i,nsrf) .LT. epsfra) trmb1(i,nsrf)=s_trmb1(i)
1844	IF (pctsrf(i,nsrf) .LT. epsfra) trmb2(i,nsrf)=s_trmb2(i)
1845	IF (pctsrf(i,nsrf) .LT. epsfra) trmb3(i,nsrf)=s_trmb3(i)
1846	ENDDO
1847	ENDDO
1848	c
1849	c
1850	c Calculer la derive du flux infrarouge
1851	c
1852	cXXX DO nsrf = 1, nbsrf
1853	DO i = 1, klon
1854	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1855	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1856	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1857	cXXX . *pctsrf(i,nsrf)
1858	cXXX ENDIF
1859	cXXX ENDDO
1860	ENDDO
1861	c
1862	c Appeler la convection (au choix)
1863	c
1864	DO k = 1, klev
1865	DO i = 1, klon
1866	conv_q(i,k) = d_q_dyn(i,k)
1867	. + d_q_vdf(i,k)/dtime
1868	conv_t(i,k) = d_t_dyn(i,k)
1869	. + d_t_vdf(i,k)/dtime
1870	ENDDO
1871	ENDDO
1872	IF (check) THEN
1873	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1874	WRITE(lunout,*) "avantcon=", za
1875	ENDIF
1876	zx_ajustq = .FALSE.
1877	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1878	IF (zx_ajustq) THEN
1879	DO i = 1, klon
1880	z_avant(i) = 0.0
1881	ENDDO
1882	DO k = 1, klev
1883	DO i = 1, klon
1884	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1885	. *(paprs(i,k)-paprs(i,k+1))/RG
1886	ENDDO
1887	ENDDO
1888	ENDIF
1889	IF (iflag_con.EQ.1) THEN
1890	stop'reactiver le call conlmd dans physiq.F'
1891	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1892	c . d_t_con, d_q_con,
1893	c . rain_con, snow_con, ibas_con, itop_con)
1894	ELSE IF (iflag_con.EQ.2) THEN
1895	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1896	e conv_t, conv_q, zxfluxq(1,1), omega,
1897	s d_t_con, d_q_con, rain_con, snow_con,
1898	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1899	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1900	WHERE (rain_con < 0.) rain_con = 0.
1901	WHERE (snow_con < 0.) snow_con = 0.
1902	DO i = 1, klon
1903	ibas_con(i) = klev+1 - kcbot(i)
1904	itop_con(i) = klev+1 - kctop(i)
1905	ENDDO
1906	ELSE IF (iflag_con.GE.3) THEN
1907	c nb of tracers for the KE convection:
1908	c MAF la partie traceurs est faite dans phytrac
1909	c on met ntra=1 pour limiter les appels mais on peut
1910	c supprimer les calculs / ftra.
1911	ntra = 1
1912	c sb, oct02:
1913	c Schema de convection modularise et vectorise:
1914	c (driver commun aux versions 3 et 4)
1915	c
1916	IF (ok_cvl) THEN ! new driver for convectL
1917
1918	CALL concvl (iflag_con,
1919	. dtime,paprs,pplay,t_seri,q_seri,
1920	. u_seri,v_seri,tr_seri,ntra,
1921	. ema_work1,ema_work2,
1922	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1923	. rain_con, snow_con, ibas_con, itop_con,
1924	. upwd,dnwd,dnwd0,
1925	. Ma,cape,tvp,iflagctrl,
1926	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
1927	. pmflxr,pmflxs,
1928	. da,phi,mp)
1929
1930	cIM cf. FH
1931	clwcon0=qcondc
1932	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
1933
1934	ELSE ! ok_cvl
1935	c MAF conema3 ne contient pas les traceurs
1936	CALL conema3 (dtime,
1937	. paprs,pplay,t_seri,q_seri,
1938	. u_seri,v_seri,tr_seri,ntra,
1939	. ema_work1,ema_work2,
1940	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1941	. rain_con, snow_con, ibas_con, itop_con,
1942	. upwd,dnwd,dnwd0,bas,top,
1943	. Ma,cape,tvp,rflag,
1944	. pbase
1945	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1946	. ,clwcon0)
1947
1948	ENDIF ! ok_cvl
1949
1950	IF (.NOT. ok_gust) THEN
1951	do i = 1, klon
1952	wd(i)=0.0
1953	enddo
1954	ENDIF
1955
1956	c =================================================================== c
1957	c Calcul des proprietes des nuages convectifs
1958	c
1959	DO k = 1, klev
1960	DO i = 1, klon
1961	zx_t = t_seri(i,k)
1962	IF (thermcep) THEN
1963	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1964	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1965	zx_qs = MIN(0.5,zx_qs)
1966	zcor = 1./(1.-retv*zx_qs)
1967	zx_qs = zx_qs*zcor
1968	ELSE
1969	IF (zx_t.LT.t_coup) THEN
1970	zx_qs = qsats(zx_t)/pplay(i,k)
1971	ELSE
1972	zx_qs = qsatl(zx_t)/pplay(i,k)
1973	ENDIF
1974	ENDIF
1975	zqsat(i,k)=zx_qs
1976	ENDDO
1977	ENDDO
1978
1979	c calcul des proprietes des nuages convectifs
1980	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1981	call clouds_gno
1982	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1983
1984	c =================================================================== c
1985
1986	DO i = 1, klon
1987	ema_pcb(i) = pbase(i)
1988	ENDDO
1989	DO i = 1, klon
1990	ema_pct(i) = paprs(i,itop_con(i))
1991	ENDDO
1992	DO i = 1, klon
1993	ema_cbmf(i) = ema_workcbmf(i)
1994	ENDDO
1995	ELSE
1996	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
1997	CALL abort
1998	ENDIF
1999
2000	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2001	c . d_u_con, d_v_con)
2002
2003	DO k = 1, klev
2004	DO i = 1, klon
2005	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
2006	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
2007	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
2008	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
2009	ENDDO
2010	ENDDO
2011	c
2012	IF (if_ebil.ge.2) THEN
2013	ztit='after convect'
2014	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2015	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2016	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2017	call diagphy(airephy,ztit,ip_ebil
2018	e , zero_v, zero_v, zero_v, zero_v, zero_v
2019	e , zero_v, rain_con, snow_con, ztsol
2020	e , d_h_vcol, d_qt, d_ec
2021	s , fs_bound, fq_bound )
2022	END IF
2023	C
2024	IF (check) THEN
2025	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2026	WRITE(lunout,*)"aprescon=", za
2027	zx_t = 0.0
2028	za = 0.0
2029	DO i = 1, klon
2030	za = za + airephy(i)/FLOAT(klon)
2031	zx_t = zx_t + (rain_con(i)+
2032	. snow_con(i))*airephy(i)/FLOAT(klon)
2033	ENDDO
2034	zx_t = zx_t/za*dtime
2035	WRITE(lunout,*)"Precip=", zx_t
2036	ENDIF
2037	IF (zx_ajustq) THEN
2038	DO i = 1, klon
2039	z_apres(i) = 0.0
2040	ENDDO
2041	DO k = 1, klev
2042	DO i = 1, klon
2043	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2044	. *(paprs(i,k)-paprs(i,k+1))/RG
2045	ENDDO
2046	ENDDO
2047	DO i = 1, klon
2048	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2049	. /z_apres(i)
2050	ENDDO
2051	DO k = 1, klev
2052	DO i = 1, klon
2053	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2054	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2055	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2056	ENDIF
2057	ENDDO
2058	ENDDO
2059	ENDIF
2060	zx_ajustq=.FALSE.
2061	c
2062	c===================================================================
2063	c Convection seche (thermiques ou ajustement)
2064	c===================================================================
2065	c
2066	d_t_ajs(:,:)=0.
2067	d_u_ajs(:,:)=0.
2068	d_v_ajs(:,:)=0.
2069	d_q_ajs(:,:)=0.
2070	fm_therm(:,:)=0.
2071	entr_therm(:,:)=0.
2072	c
2073	IF(prt_level>9)WRITE(lunout,*)
2074	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2075	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2076	if(iflag_thermals.lt.0) then
2077	c Rien
2078	c ====
2079	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2080	else if(iflag_thermals.eq.0) then
2081
2082	c Ajustement sec
2083	c ==============
2084	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2085	CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
2086	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
2087	q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
2088	else
2089	c Thermiques
2090	c ==========
2091	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2092	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2093	call calltherm(pdtphys
2094	s ,pplay,paprs,pphi
2095	s ,u_seri,v_seri,t_seri,q_seri
2096	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2097	s ,fm_therm,entr_therm)
2098	endif
2099	c
2100	c===================================================================
2101	c
2102	IF (if_ebil.ge.2) THEN
2103	ztit='after dry_adjust'
2104	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2105	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2106	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2107	END IF
2108
2109
2110	c-------------------------------------------------------------------------
2111	c Caclul des ratqs
2112	c-------------------------------------------------------------------------
2113
2114	c print*,'calcul des ratqs'
2115	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2116	c ----------------
2117	c on ecrase le tableau ratqsc calcule par clouds_gno
2118	if (iflag_cldcon.eq.1) then
2119	do k=1,klev
2120	do i=1,klon
2121	if(ptconv(i,k)) then
2122	ratqsc(i,k)=ratqsbas
2123	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2124	else
2125	ratqsc(i,k)=0.
2126	endif
2127	enddo
2128	enddo
2129	endif
2130
2131	c ratqs stables
2132	c -------------
2133	do k=1,klev
2134	cIM RAJOUT boucle do=i
2135	do i=1, klon
2136	cIM ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
2137	cIM s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
2138	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2139	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2140	cIM print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
2141	enddo
2142	enddo
2143
2144
2145	c ratqs final
2146	c -----------
2147	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
2148	c les ratqs sont une conbinaison de ratqss et ratqsc
2149	c ratqs final
2150	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
2151	c relaxation des ratqs
2152	c facttemps=exp(-pdtphys/1.e4)
2153	facteur=exp(-pdtphys*facttemps)
2154	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
2155	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
2156	c print*,'calcul des ratqs fini'
2157	else
2158	c on ne prend que le ratqs stable pour fisrtilp
2159	ratqs(:,:)=ratqss(:,:)
2160	endif
2161
2162
2163	c
2164	c Appeler le processus de condensation a grande echelle
2165	c et le processus de precipitation
2166	c-------------------------------------------------------------------------
2167	CALL fisrtilp(dtime,paprs,pplay,
2168	. t_seri, q_seri,ptconv,ratqs,
2169	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2170	. rain_lsc, snow_lsc,
2171	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2172	. frac_impa, frac_nucl,
2173	. prfl, psfl, rhcl)
2174
2175	WHERE (rain_lsc < 0) rain_lsc = 0.
2176	WHERE (snow_lsc < 0) snow_lsc = 0.
2177	DO k = 1, klev
2178	DO i = 1, klon
2179	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2180	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2181	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2182	cldfra(i,k) = rneb(i,k)
2183	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2184	ENDDO
2185	ENDDO
2186	IF (check) THEN
2187	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2188	WRITE(lunout,*)"apresilp=", za
2189	zx_t = 0.0
2190	za = 0.0
2191	DO i = 1, klon
2192	za = za + airephy(i)/FLOAT(klon)
2193	zx_t = zx_t + (rain_lsc(i)
2194	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2195	ENDDO
2196	zx_t = zx_t/za*dtime
2197	WRITE(lunout,*)"Precip=", zx_t
2198	ENDIF
2199	c
2200	IF (if_ebil.ge.2) THEN
2201	ztit='after fisrt'
2202	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2203	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2204	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2205	call diagphy(airephy,ztit,ip_ebil
2206	e , zero_v, zero_v, zero_v, zero_v, zero_v
2207	e , zero_v, rain_lsc, snow_lsc, ztsol
2208	e , d_h_vcol, d_qt, d_ec
2209	s , fs_bound, fq_bound )
2210	END IF
2211	c
2212	c-------------------------------------------------------------------
2213	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2214	c-------------------------------------------------------------------
2215
2216	c 1. NUAGES CONVECTIFS
2217	c
2218	cIM cf FH
2219	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2220	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2221	snow_tiedtke=0.
2222	c print*,'avant calcul de la pseudo precip '
2223	c print*,'iflag_cldcon',iflag_cldcon
2224	if (iflag_cldcon.eq.-1) then
2225	rain_tiedtke=rain_con
2226	else
2227	c print*,'calcul de la pseudo precip '
2228	rain_tiedtke=0.
2229	c print*,'calcul de la pseudo precip 0'
2230	do k=1,klev
2231	do i=1,klon
2232	if (d_q_con(i,k).lt.0.) then
2233	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2234	s *(paprs(i,k)-paprs(i,k+1))/rg
2235	endif
2236	enddo
2237	enddo
2238	endif
2239	c
2240	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2241	c
2242
2243	c Nuages diagnostiques pour Tiedtke
2244	CALL diagcld1(paprs,pplay,
2245	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2246	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2247	. diafra,dialiq)
2248	DO k = 1, klev
2249	DO i = 1, klon
2250	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2251	cldliq(i,k) = dialiq(i,k)
2252	cldfra(i,k) = diafra(i,k)
2253	ENDIF
2254	ENDDO
2255	ENDDO
2256
2257	ELSE IF (iflag_cldcon.eq.3) THEN
2258	c On prend pour les nuages convectifs le max du calcul de la
2259	c convection et du calcul du pas de temps précédent diminué d'un facteur
2260	c facttemps
2261	c facttemps=pdtphys/1.e4
2262	facteur = pdtphys *facttemps
2263	do k=1,klev
2264	do i=1,klon
2265	rnebcon(i,k)=rnebcon(i,k)*facteur
2266	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2267	s then
2268	rnebcon(i,k)=rnebcon0(i,k)
2269	clwcon(i,k)=clwcon0(i,k)
2270	endif
2271	enddo
2272	enddo
2273
2274	c
2275	cIM calcul nuages par le simulateur ISCCP
2276	c
2277	IF (ok_isccp) THEN
2278	#include "calcul_simulISCCP.h"
2279	ENDIF !ok_isccp
2280
2281	c On prend la somme des fractions nuageuses et des contenus en eau
2282	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2283	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2284
2285	ENDIF
2286
2287	c
2288	c 2. NUAGES STARTIFORMES
2289	c
2290	IF (ok_stratus) THEN
2291	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2292	DO k = 1, klev
2293	DO i = 1, klon
2294	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2295	cldliq(i,k) = dialiq(i,k)
2296	cldfra(i,k) = diafra(i,k)
2297	ENDIF
2298	ENDDO
2299	ENDDO
2300	ENDIF
2301	c
2302	c Precipitation totale
2303	c
2304	DO i = 1, klon
2305	rain_fall(i) = rain_con(i) + rain_lsc(i)
2306	snow_fall(i) = snow_con(i) + snow_lsc(i)
2307	ENDDO
2308	c
2309	IF (if_ebil.ge.2) THEN
2310	ztit="after diagcld"
2311	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2312	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2313	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2314	END IF
2315	c
2316	c Calculer l'humidite relative pour diagnostique
2317	c
2318	DO k = 1, klev
2319	DO i = 1, klon
2320	zx_t = t_seri(i,k)
2321	IF (thermcep) THEN
2322	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2323	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2324	zx_qs = MIN(0.5,zx_qs)
2325	zcor = 1./(1.-retv*zx_qs)
2326	zx_qs = zx_qs*zcor
2327	ELSE
2328	IF (zx_t.LT.t_coup) THEN
2329	zx_qs = qsats(zx_t)/pplay(i,k)
2330	ELSE
2331	zx_qs = qsatl(zx_t)/pplay(i,k)
2332	ENDIF
2333	ENDIF
2334	zx_rh(i,k) = q_seri(i,k)/zx_qs
2335	zqsat(i,k)=zx_qs
2336	ENDDO
2337	ENDDO
2338	cjq - introduce the aerosol direct and first indirect radiative forcings
2339	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2340	IF (ok_ade.OR.ok_aie) THEN
2341	! Get sulfate aerosol distribution
2342	CALL readsulfate(rjourvrai, debut, sulfate)
2343	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2344
2345	! Calculate aerosol optical properties (Olivier Boucher)
2346	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2347	. tau_ae, piz_ae, cg_ae, aerindex)
2348	cym
2349	ELSE
2350	tau_ae(:,:,:)=0.0
2351	piz_ae(:,:,:)=0.0
2352	cg_ae(:,:,:)=0.0
2353	cym
2354	ENDIF
2355
2356	#ifdef INCA
2357	calday = FLOAT(julien) + gmtime
2358
2359	#ifdef INCA_AER
2360	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2361	& prfl,psfl,pctsrf(1,3),airephy,xjour,rlat,rlon)
2362	#endif
2363
2364	#ifdef INCAINFO
2365	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2366	#endif
2367
2368	CALL chemhook_begin (calday,
2369	$ pctsrf(1,1),
2370	$ rlat,
2371	$ rlon,
2372	$ airephy,
2373	$ paprs,
2374	$ pplay,
2375	$ ycoefh,
2376	$ pphi,
2377	$ t_seri,
2378	$ u,
2379	$ v,
2380	$ wo,
2381	$ q_seri,
2382	$ zxtsol,
2383	$ zxsnow,
2384	$ solsw,
2385	$ albsol,
2386	$ rain_fall,
2387	$ snow_fall,
2388	$ itop_con,
2389	$ ibas_con,
2390	$ cldfra,
2391	$ iim,
2392	$ jjm,
2393	#ifdef INCA_AER
2394	$ tr_seri,
2395	$ ftsol,
2396	$ paprs,
2397	$ cdragh,
2398	$ cdragm,
2399	$ pctsrf,
2400	$ pdtphys,
2401	$ itap)
2402	#else
2403	$ tr_seri)
2404	#endif
2405
2406
2407	#ifdef INCAINFO
2408	WRITE(lunout,*)'OK.'
2409	#endif
2410	#endif
2411	c
2412	c Calculer les parametres optiques des nuages et quelques
2413	c parametres pour diagnostiques:
2414	c
2415	if (ok_newmicro) then
2416	CALL newmicro (paprs, pplay,ok_newmicro,
2417	. t_seri, cldliq, cldfra, cldtau, cldemi,
2418	. cldh, cldl, cldm, cldt, cldq,
2419	. flwp, fiwp, flwc, fiwc,
2420	e ok_aie,
2421	e sulfate, sulfate_pi,
2422	e bl95_b0, bl95_b1,
2423	s cldtaupi, re, fl)
2424	else
2425	CALL nuage (paprs, pplay,
2426	. t_seri, cldliq, cldfra, cldtau, cldemi,
2427	. cldh, cldl, cldm, cldt, cldq,
2428	e ok_aie,
2429	e sulfate, sulfate_pi,
2430	e bl95_b0, bl95_b1,
2431	s cldtaupi, re, fl)
2432
2433	endif
2434	c
2435	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2436	c
2437	IF (MOD(itaprad,radpas).EQ.0) THEN
2438	DO i = 1, klon
2439	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2440	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2441	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2442	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2443	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2444	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2445	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2446	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2447	ENDDO
2448	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2449	e (dist, rmu0, fract,
2450	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2451	e wo,
2452	e cldfra, cldemi, cldtau,
2453	s heat,heat0,cool,cool0,radsol,albpla,
2454	s topsw,toplw,solsw,sollw,
2455	s sollwdown,
2456	s topsw0,toplw0,solsw0,sollw0,
2457	s lwdn0, lwdn, lwup0, lwup,
2458	s swdn0, swdn, swup0, swup,
2459	e ok_ade, ok_aie, ! new for aerosol radiative effects
2460	e tau_ae, piz_ae, cg_ae, ! ="=
2461	s topswad, solswad, ! ="=
2462	e cldtaupi, ! ="=
2463	s topswai, solswai) ! ="=
2464	itaprad = 0
2465	ENDIF
2466	itaprad = itaprad + 1
2467	c
2468	c Ajouter la tendance des rayonnements (tous les pas)
2469	c
2470	DO k = 1, klev
2471	DO i = 1, klon
2472	t_seri(i,k) = t_seri(i,k)
2473	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2474	ENDDO
2475	ENDDO
2476	c
2477	IF (if_ebil.ge.2) THEN
2478	ztit='after rad'
2479	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2480	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2481	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2482	call diagphy(airephy,ztit,ip_ebil
2483	e , topsw, toplw, solsw, sollw, zero_v
2484	e , zero_v, zero_v, zero_v, ztsol
2485	e , d_h_vcol, d_qt, d_ec
2486	s , fs_bound, fq_bound )
2487	END IF
2488	c
2489	c
2490	c Calculer l'hydrologie de la surface
2491	c
2492	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2493	c . agesno, ftsol,fqsurf,fsnow, ruis)
2494	c
2495	DO i = 1, klon
2496	zxqsurf(i) = 0.0
2497	zxsnow(i) = 0.0
2498	ENDDO
2499	DO nsrf = 1, nbsrf
2500	DO i = 1, klon
2501	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2502	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2503	ENDDO
2504	ENDDO
2505	c
2506	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2507	c
2508	cXXX DO nsrf = 1, nbsrf
2509	cXXX DO i = 1, klon
2510	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
2511	cXXX fqsurf(i,nsrf) = zxqsurf(i)
2512	cXXX fsnow(i,nsrf) = zxsnow(i)
2513	cXXX ENDIF
2514	cXXX ENDDO
2515	cXXX ENDDO
2516	c
2517	c Calculer le bilan du sol et la derive de temperature (couplage)
2518	c
2519	DO i = 1, klon
2520	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2521	c a la demande de JLD
2522	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2523	ENDDO
2524	c
2525	cmoddeblott(jan95)
2526	c Appeler le programme de parametrisation de l'orographie
2527	c a l'echelle sous-maille:
2528	c
2529	IF (ok_orodr) THEN
2530	c
2531	c selection des points pour lesquels le shema est actif:
2532	igwd=0
2533	DO i=1,klon
2534	itest(i)=0
2535	c IF ((zstd(i).gt.10.0)) THEN
2536	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2537	itest(i)=1
2538	igwd=igwd+1
2539	idx(igwd)=i
2540	ENDIF
2541	ENDDO
2542	c igwdim=MAX(1,igwd)
2543	c
2544	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2545	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2546	e igwd,idx,itest,
2547	e t_seri, u_seri, v_seri,
2548	cIM 141004 s zulow, zvlow, zustr, zvstr,
2549	s zulow, zvlow, zustrdr, zvstrdr,
2550	s d_t_oro, d_u_oro, d_v_oro)
2551	c
2552	c ajout des tendances
2553	DO k = 1, klev
2554	DO i = 1, klon
2555	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2556	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2557	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2558	ENDDO
2559	ENDDO
2560	c
2561	ENDIF ! fin de test sur ok_orodr
2562	c
2563	IF (ok_orolf) THEN
2564	c
2565	c selection des points pour lesquels le shema est actif:
2566	igwd=0
2567	DO i=1,klon
2568	itest(i)=0
2569	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2570	itest(i)=1
2571	igwd=igwd+1
2572	idx(igwd)=i
2573	ENDIF
2574	ENDDO
2575	c igwdim=MAX(1,igwd)
2576	c
2577	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2578	e rlat,zmea,zstd,zpic,
2579	e itest,
2580	e t_seri, u_seri, v_seri,
2581	s zulow, zvlow, zustrli, zvstrli,
2582	s d_t_lif, d_u_lif, d_v_lif)
2583	c
2584	c ajout des tendances
2585	DO k = 1, klev
2586	DO i = 1, klon
2587	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2588	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2589	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2590	ENDDO
2591	ENDDO
2592	c
2593	ENDIF ! fin de test sur ok_orolf
2594	c
2595	cIM cf. FLott BEG
2596	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
2597
2598	DO i = 1, klon
2599	zustrph(i)=0.
2600	zvstrph(i)=0.
2601	ENDDO
2602	DO k = 1, klev
2603	DO i = 1, klon
2604	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
2605	c (paprs(i,k)-paprs(i,k+1))/rg
2606	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
2607	c (paprs(i,k)-paprs(i,k+1))/rg
2608	ENDDO
2609	ENDDO
2610	c
2611	cIM calcul composantes axiales du moment angulaire et couple des montagnes
2612	c
2613	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
2614	C ra,rg,romega,
2615	C rlat,rlon,pphis,
2616	C zustrdr,zustrli,zustrph,
2617	C zvstrdr,zvstrli,zvstrph,
2618	C paprs,u,v,
2619	C aam, torsfc)
2620	cIM cf. FLott END
2621	c
2622	IF (if_ebil.ge.2) THEN
2623	ztit='after orography'
2624	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2625	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2626	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2627	END IF
2628	c
2629	c
2630	cAA
2631	cAA Installation de l'interface online-offline pour traceurs
2632	cAA
2633	c====================================================================
2634	c Calcul des tendances traceurs
2635	c====================================================================
2636	C
2637	call phytrac (rnpb,
2638	I itap,
2639	I julien,
2640	I gmtime,
2641	I debut,
2642	I lafin,
2643	I nqmax-2,
2644	I nlon,
2645	I nlev,
2646	I dtime,
2647	I u,
2648	I v,
2649	I t,
2650	I paprs,
2651	I pplay,
2652	I pmfu,
2653	I pmfd,
2654	I pen_u,
2655	I pde_u,
2656	I pen_d,
2657	I pde_d,
2658	I ycoefh,
2659	I fm_therm,
2660	I entr_therm,
2661	I yu1,
2662	I yv1,
2663	I ftsol,
2664	I pctsrf,
2665	I rlat,
2666	I frac_impa,
2667	I frac_nucl,
2668	I rlon,
2669	I presnivs,
2670	I pphis,
2671	I pphi,
2672	I albsol,
2673	I qx(1,1,1),
2674	I rhcl,
2675	I cldfra,
2676	I rneb,
2677	I diafra,
2678	I cldliq,
2679	I itop_con,
2680	I ibas_con,
2681	I pmflxr,
2682	I pmflxs,
2683	I prfl,
2684	I psfl,
2685	I da,
2686	I phi,
2687	I mp,
2688	I upwd,
2689	I dnwd,
2690	#ifdef INCA
2691	I flxmass_w,
2692	#endif
2693	O tr_seri)
2694
2695	IF (offline) THEN
2696
2697	print*,'Attention on met a 0 les thermiques pour phystoke'
2698	call phystokenc (
2699	I nlon,nlev,pdtphys,rlon,rlat,
2700	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2701	I fm_therm,entr_therm,
2702	I ycoefh,yu1,yv1,ftsol,pctsrf,
2703	I frac_impa, frac_nucl,
2704	I pphis,airephy,dtime,itap)
2705
2706
2707	ENDIF
2708
2709	c
2710	c Calculer le transport de l'eau et de l'energie (diagnostique)
2711	c
2712	CALL transp (paprs,zxtsol,
2713	e t_seri, q_seri, u_seri, v_seri, zphi,
2714	s ve, vq, ue, uq)
2715	c
2716	cIM diag. bilKP
2717	c
2718	CALL transp_lay (paprs,zxtsol,
2719	e t_seri, q_seri, u_seri, v_seri, zphi,
2720	s ve_lay, vq_lay, ue_lay, uq_lay)
2721	c
2722	c Accumuler les variables a stocker dans les fichiers histoire:
2723	c
2724	c+jld ec_conser
2725	DO k = 1, klev
2726	DO i = 1, klon
2727	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
2728	d_t_ec(i,k)=0.5/ZRCPD
2729	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
2730	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
2731	d_t_ec(i,k) = d_t_ec(i,k)/dtime
2732	END DO
2733	END DO
2734	c-jld ec_conser
2735	IF (if_ebil.ge.1) THEN
2736	ztit='after physic'
2737	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
2738	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2739	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2740	C Comme les tendances de la physique sont ajoute dans la dynamique,
2741	C on devrait avoir que la variation d'entalpie par la dynamique
2742	C est egale a la variation de la physique au pas de temps precedent.
2743	C Donc la somme de ces 2 variations devrait etre nulle.
2744	call diagphy(airephy,ztit,ip_ebil
2745	e , topsw, toplw, solsw, sollw, sens
2746	e , evap, rain_fall, snow_fall, ztsol
2747	e , d_h_vcol, d_qt, d_ec
2748	s , fs_bound, fq_bound )
2749	C
2750	d_h_vcol_phy=d_h_vcol
2751	C
2752	END IF
2753	C
2754	c=======================================================================
2755	c SORTIES
2756	c=======================================================================
2757
2758	cIM Interpolation sur les niveaux de pression du NMC
2759	c -------------------------------------------------
2760	c
2761	#include "calcul_STDlev.h"
2762	c
2763	c slp sea level pressure
2764	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
2765	c
2766	ccc prw = eau precipitable
2767	DO i = 1, klon
2768	prw(i) = 0.
2769	DO k = 1, klev
2770	prw(i) = prw(i) +
2771	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
2772	ENDDO
2773	ENDDO
2774	c
2775	cIM initialisation + calculs divers diag AMIP2
2776	c
2777	#include "calcul_divers.h"
2778	c
2779	#ifdef INCA
2780	#ifdef INCAINFO
2781	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
2782	#endif
2783	CALL chemhook_end (calday,
2784	$ dtime,
2785	$ pplay,
2786	$ t_seri,
2787	$ tr_seri,
2788	$ nbtr,
2789	$ paprs,
2790	#ifdef INCA_CH4
2791	$ q_seri,
2792	#endif
2793	$ annee_ref,
2794	$ day_ini,
2795	#ifdef INCA_AER
2796	$ xjour,
2797	$ pphi,
2798	$ pphis,
2799	$ zx_rh,
2800	$ qx(1,1,1))
2801	#else
2802	$ xjour)
2803	#endif
2804	#ifdef INCAINFO
2805	WRITE(lunout,*)'OK.'
2806	#endif
2807	#endif
2808
2809	c=============================================================
2810	c
2811	c Convertir les incrementations en tendances
2812	c
2813	DO k = 1, klev
2814	DO i = 1, klon
2815	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2816	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2817	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2818	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2819	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2820	ENDDO
2821	ENDDO
2822	c
2823	IF (nqmax.GE.3) THEN
2824	DO iq = 3, nqmax
2825	DO k = 1, klev
2826	DO i = 1, klon
2827	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2828	ENDDO
2829	ENDDO
2830	ENDDO
2831	ENDIF
2832	c
2833	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
2834	c#include "write_bilKP_ins.h"
2835	c#include "write_bilKP_ave.h"
2836	c
2837	c Sauvegarder les valeurs de t et q a la fin de la physique:
2838	c
2839	DO k = 1, klev
2840	DO i = 1, klon
2841	t_ancien(i,k) = t_seri(i,k)
2842	q_ancien(i,k) = q_seri(i,k)
2843	ENDDO
2844	ENDDO
2845	c
2846	c 22.03.04 BEG
2847	c=============================================================
2848	c Ecriture des sorties
2849	c=============================================================
2850	#ifdef CPP_IOIPSL
2851
2852	#ifdef histhf
2853	#include "write_histhf.h"
2854	#endif
2855
2856	#ifdef histday
2857	#include "write_histday.h"
2858	#include "write_histday_seri.h"
2859	#endif
2860
2861	#ifdef histmth
2862	#include "write_histmth.h"
2863	#endif
2864
2865	#ifdef histins
2866	#include "write_histins.h"
2867	#endif
2868
2869	#ifdef histISCCP
2870	#include "write_histISCCP.h"
2871	#endif
2872
2873	#ifdef histmthNMC
2874	#include "write_histmthNMC.h"
2875	#endif
2876
2877	#endif
2878
2879	c 22.03.04 END
2880	c
2881	c====================================================================
2882	c Si c'est la fin, il faut conserver l'etat de redemarrage
2883	c====================================================================
2884	c
2885	IF (lafin) THEN
2886	itau_phy = itau_phy + itap
2887	ccc IF (ok_oasis) CALL quitcpl
2888	CALL phyredem ("restartphy.nc",dtime,radpas,
2889	. rlat, rlon, pctsrf, ftsol, ftsoil,
2890	cIM "slab" ocean
2891	. tslab, seaice,
2892	. fqsurf, qsol,
2893	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
2894	. solsw, sollwdown,dlw,
2895	. radsol,frugs,agesno,
2896	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2897	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
2898	ENDIF
2899
2900
2901	RETURN
2902	END
2903	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2904	IMPLICIT none
2905	c
2906	c Calculer et imprimer l'eau totale. A utiliser pour verifier
2907	c la conservation de l'eau
2908	c
2909	#include "YOMCST.h"
2910	INTEGER klon,klev
2911	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2912	REAL aire(klon)
2913	REAL qtotal, zx, qcheck
2914	INTEGER i, k
2915	c
2916	zx = 0.0
2917	DO i = 1, klon
2918	zx = zx + aire(i)
2919	ENDDO
2920	qtotal = 0.0
2921	DO k = 1, klev
2922	DO i = 1, klon
2923	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2924	. *(paprs(i,k)-paprs(i,k+1))/RG
2925	ENDDO
2926	ENDDO
2927	c
2928	qcheck = qtotal/zx
2929	c
2930	RETURN
2931	END
2932	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2933	IMPLICIT none
2934	c
2935	c Tranformer une variable de la grille physique a
2936	c la grille d'ecriture
2937	c
2938	INTEGER nfield,nlon,iim,jjmp1, jjm
2939	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2940	c
2941	INTEGER i, n, ig
2942	c
2943	jjm = jjmp1 - 1
2944	DO n = 1, nfield
2945	DO i=1,iim
2946	ecrit(i,n) = fi(1,n)
2947	ecrit(i+jjm*iim,n) = fi(nlon,n)
2948	ENDDO
2949	DO ig = 1, nlon - 2
2950	ecrit(iim+ig,n) = fi(1+ig,n)
2951	ENDDO
2952	ENDDO
2953	RETURN
2954	END

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

Download in other formats:

Original Format