Context Navigation

physiq.F @ 5373

Last change on this file since 5373 was 755, checked in by lsce, 18 years ago
Ajout une variable d'entree pour une routine inca ACo
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 117.5 KB

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

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

Download in other formats:

Original Format