Context Navigation

physiq.F @ 1643

Last change on this file since 1643 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

physiq.F : appel modifie a calltherm et fisrtilp
thermcell_main.F90 : changement dans les flag_thermals_ed
thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

cloudth.F90 : PDF d'eau bi-gaussiennes

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

File size: 123.6 KB

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

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

Download in other formats:

Original Format