Context Navigation

physiq.F @ 979

Last change on this file since 979 was 979, checked in by Laurent Fairhead, 16 years ago
Inclusion d'un seuil de validite sur ape et ale pour les poches froides FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 110.1 KB

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

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

Download in other formats:

Original Format