Context Navigation

physiq.F @ 839

Last change on this file since 839 was 839, checked in by lmdzadmin, 17 years ago
Appel simulateur ISCCP au cas ou "#define histISCCP" IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 107.0 KB

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

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

Download in other formats:

Original Format