Context Navigation

physiq.F @ 1287

Last change on this file since 1287 was 954, checked in by lsce, 17 years ago

ACo+JG

Ajout du flag aerosol_couple :
false=lecture des sulfates dans un fichier(par defaut) - configuration existante
true=calcul des aerosol par INCA

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

File size: 108.1 KB

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

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

Download in other formats:

Original Format