Context Navigation

physiq.F @ 1046

Last change on this file since 1046 was 1045, checked in by lmdzadmin, 16 years ago
On sort meantaucld dans histISCCP, car sortie simulateur (enleve du histday) Lecture fichiers input ISCCP "debut" physiq IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 113.9 KB

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

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

Download in other formats:

Original Format