Context Navigation

physiq.F @ 1059

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

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

Download in other formats:

Original Format