Context Navigation

physiq.F @ 1533

Last change on this file since 1533 was 1533, checked in by musat, 13 years ago

Ajout fichier 1D mensuel paramLMDZ_phy.nc contenant

les parametres orbitaux
les taux des GES
les moyennes globales ((bils, evap, evap_land, flat, nettop0, nettop, precip, tsol, t2m, prw)

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

File size: 120.6 KB

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

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

Download in other formats:

Original Format