Context Navigation

physiq.F @ 1373

Last change on this file since 1373 was 1373, checked in by musat, 14 years ago

Some changes for new physics :

new formulation for wakes' alp_offset (physiq)
bug fix for lalim in thermals plume
add iflag_pbl=9 option (yamada4)

FH/IM

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

File size: 123.4 KB

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

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

Download in other formats:

Original Format