Context Navigation

physiq.F @ 717

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

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

Download in other formats:

Original Format