Context Navigation

physiq.F @ 5172

Last change on this file since 5172 was 1428, checked in by Laurent Fairhead, 14 years ago

Problems with some initialisations in the new physics package. These are
temporary fixes. Further investigations needed. Trac ticket issued. JYG

Problèmes d'initialisations dans la nouvelle physique. Ces corrections
sont temporaires, il faut trouver la vraie raison des plantages. Un ticket
trac est créé. JYG

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

File size: 123.0 KB

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

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

Download in other formats:

Original Format