Context Navigation

physiq.F @ 1457

Last change on this file since 1457 was 1457, checked in by musat, 14 years ago
Ooops phyetat0, phyredem: correction dimension verticale pbl_tke: pbl_tke(:,1:klev+1,:) physiq: pour pouvoir fixer la longitude solaire avec la nouvelle orbite
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 123.1 KB

Line
1	! $Id: physiq.F 1457 2010-11-22 15:19:25Z musat $
2	c#define IO_DEBUG
3
4	SUBROUTINE physiq (nlon,nlev,
5	. debut,lafin,jD_cur, jH_cur,pdtphys,
6	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
7	. u,v,t,qx,
8	. flxmass_w,
9	. d_u, d_v, d_t, d_qx, d_ps
10	. , dudyn
11	. , PVteta)
12
13	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
14	$ histwrite, ju2ymds, ymds2ju, ioget_year_len
15	USE comgeomphy
16	USE phys_cal_mod
17	USE write_field_phy
18	USE dimphy
19	USE infotrac
20	USE mod_grid_phy_lmdz
21	USE mod_phys_lmdz_para
22	USE iophy
23	USE misc_mod, mydebug=>debug
24	USE vampir
25	USE pbl_surface_mod, ONLY : pbl_surface
26	USE change_srf_frac_mod
27	USE surface_data, ONLY : type_ocean, ok_veget
28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
29	USE phys_state_var_mod ! Variables sauvegardees de la physique
30	USE phys_output_var_mod ! Variables pour les ecritures des sorties
31	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
32	USE phys_output_mod
33	use open_climoz_m, only: open_climoz ! ozone climatology from a file
34	use regr_pr_av_m, only: regr_pr_av
35	use netcdf95, only: nf95_close
36	cIM for NMC files
37	use netcdf, only: nf90_fill_real
38	use mod_phys_lmdz_mpi_data, only: is_mpi_root
39	USE aero_mod
40	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
41	use conf_phys_m, only: conf_phys
42	use radlwsw_m, only: radlwsw
43	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	cIM 180608
1253	c pmflxr=0.
1254	c pmflxs=0.
1255
1256	itau_con=0
1257	first=.false.
1258
1259	endif ! first
1260
1261	modname = 'physiq'
1262	cIM
1263	IF (ip_ebil_phy.ge.1) THEN
1264	DO i=1,klon
1265	zero_v(i)=0.
1266	END DO
1267	END IF
1268	ok_sync=.TRUE.
1269
1270	IF (debut) THEN
1271	CALL suphel ! initialiser constantes et parametres phys.
1272	ENDIF
1273
1274	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1275
1276
1277	c======================================================================
1278	! Gestion calendrier : mise a jour du module phys_cal_mod
1279	!
1280	c CALL phys_cal_update(jD_cur,jH_cur)
1281
1282	c
1283	c Si c'est le debut, il faut initialiser plusieurs choses
1284	c ********
1285	c
1286	IF (debut) THEN
1287	!rv
1288	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1289	cde la convection a partir des caracteristiques du thermique
1290	wght_th(:,:)=1.
1291	lalim_conv(:)=1
1292	cRC
1293	u10m(:,:)=0.
1294	v10m(:,:)=0.
1295	rain_con(:)=0.
1296	snow_con(:)=0.
1297	topswai(:)=0.
1298	topswad(:)=0.
1299	solswai(:)=0.
1300	solswad(:)=0.
1301
1302	lambda_th(:,:)=0.
1303	wmax_th(:)=0.
1304	tau_overturning_th(:)=0.
1305
1306	IF (config_inca /= 'none') THEN
1307	! jg : initialisation jusqu'au ces variables sont dans restart
1308	ccm(:,:,:) = 0.
1309	tau_aero(:,:,:,:) = 0.
1310	piz_aero(:,:,:,:) = 0.
1311	cg_aero(:,:,:,:) = 0.
1312	END IF
1313
1314	rnebcon0(:,:) = 0.0
1315	clwcon0(:,:) = 0.0
1316	rnebcon(:,:) = 0.0
1317	clwcon(:,:) = 0.0
1318
1319	cIM
1320	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1321	c
1322	print*,'iflag_coupl,iflag_clos,iflag_wake',
1323	. iflag_coupl,iflag_clos,iflag_wake
1324	print*,'CYCLE_DIURNE', cycle_diurne
1325	c
1326	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1327	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1328	CALL abort_gcm (modname,abort_message,1)
1329	ENDIF
1330	c
1331	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1332	abort_message = 'ISCCP-like outputs may be available for KE
1333	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1334	CALL abort_gcm (modname,abort_message,1)
1335	ENDIF
1336	c
1337	c Initialiser les compteurs:
1338	c
1339	itap = 0
1340	itaprad = 0
1341
1342	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1343	!! Un petit travail à faire ici.
1344	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1345
1346	if (iflag_pbl>1) then
1347	PRINT*, "Using method MELLOR&YAMADA"
1348	endif
1349
1350	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1351	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1352	! dyn3d
1353	! Attention : la version precedente n'etait pas tres propre.
1354	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1355	! pour obtenir le meme resultat.
1356	dtime=pdtphys
1357	radpas = NINT( 86400./dtime/nbapp_rad)
1358	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1359
1360	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1361	cIM begin
1362	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1363	$,ratqs(1,1)
1364	cIM end
1365
1366
1367
1368	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1369	c
1370	C on remet le calendrier a zero
1371	c
1372	IF (raz_date .eq. 1) THEN
1373	itau_phy = 0
1374	ENDIF
1375
1376	cIM cf. AM 081204 BEG
1377	PRINT*,'cycle_diurne3 =',cycle_diurne
1378	cIM cf. AM 081204 END
1379	c
1380	CALL printflag( tabcntr0,radpas,ok_journe,
1381	, ok_instan, ok_region )
1382	c
1383	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1384	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1385	. pdtphys
1386	abort_message='Pas physique n est pas correct '
1387	! call abort_gcm(modname,abort_message,1)
1388	dtime=pdtphys
1389	ENDIF
1390	IF (nlon .NE. klon) THEN
1391	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1392	. klon
1393	abort_message='nlon et klon ne sont pas coherents'
1394	call abort_gcm(modname,abort_message,1)
1395	ENDIF
1396	IF (nlev .NE. klev) THEN
1397	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1398	. klev
1399	abort_message='nlev et klev ne sont pas coherents'
1400	call abort_gcm(modname,abort_message,1)
1401	ENDIF
1402	c
1403	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1404	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1405	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1406	abort_message='Nbre d appels au rayonnement insuffisant'
1407	call abort_gcm(modname,abort_message,1)
1408	ENDIF
1409	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1410	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1411	. ok_cvl
1412	c
1413	cKE43
1414	c Initialisation pour la convection de K.E. (sb):
1415	IF (iflag_con.GE.3) THEN
1416
1417	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1418	WRITE(lunout,*)
1419	. "On va utiliser le melange convectif des traceurs qui"
1420	WRITE(lunout,*)"est calcule dans convect4.3"
1421	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1422
1423	DO i = 1, klon
1424	ema_cbmf(i) = 0.
1425	ema_pcb(i) = 0.
1426	ema_pct(i) = 0.
1427	c ema_workcbmf(i) = 0.
1428	ENDDO
1429	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1430	DO i = 1, klon
1431	ibas_con(i) = 1
1432	itop_con(i) = 1
1433	ENDDO
1434	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1435	c===============================================================================
1436	cCR:04.12.07: initialisations poches froides
1437	c Controle de ALE et ALP pour la fermeture convective (jyg)
1438	if (iflag_wake.eq.1) then
1439	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1440	s ,alp_bl_prescr, ale_bl_prescr)
1441	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1442	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1443	endif
1444
1445	do i = 1,klon
1446	Ale_bl(i)=0.
1447	Alp_bl(i)=0.
1448	enddo
1449
1450	c================================================================================
1451
1452	ENDIF !debut
1453
1454	DO i=1,klon
1455	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1456	ENDDO
1457
1458	c34EK
1459	IF (ok_orodr) THEN
1460
1461	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1462	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1463	! justement quand ok_orodr = false.
1464	! ce rugoro est utilise par la couche limite et fait double emploi
1465	! avec les paramétrisations spécifiques de Francois Lott.
1466	! DO i=1,klon
1467	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1468	! ENDDO
1469	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1470	IF (ok_strato) THEN
1471	CALL SUGWD_strato(klon,klev,paprs,pplay)
1472	ELSE
1473	CALL SUGWD(klon,klev,paprs,pplay)
1474	ENDIF
1475
1476	DO i=1,klon
1477	zuthe(i)=0.
1478	zvthe(i)=0.
1479	if(zstd(i).gt.10.)then
1480	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1481	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1482	endif
1483	ENDDO
1484	ENDIF
1485	c
1486	c
1487	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1488	WRITE(lunout,*)'La frequence de lecture surface est de ',
1489	. lmt_pas
1490	c
1491	cIM 030306 END
1492
1493	capemaxcels = 't_max(X)'
1494	t2mincels = 't_min(X)'
1495	t2maxcels = 't_max(X)'
1496	tinst = 'inst(X)'
1497	tave = 'ave(X)'
1498	cIM cf. AM 081204 BEG
1499	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1500	cIM cf. AM 081204 END
1501	c
1502	c=============================================================
1503	c Initialisation des sorties
1504	c=============================================================
1505
1506	#ifdef CPP_IOIPSL
1507
1508	c$OMP MASTER
1509	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1510	& ctetaSTD,dtime,ok_veget,
1511	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1512	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1513	& read_climoz, new_aod, aerosol_couple
1514	& )
1515	c$OMP END MASTER
1516	c$OMP BARRIER
1517
1518	#ifdef histISCCP
1519	#include "ini_histISCCP.h"
1520	#endif
1521
1522	#ifdef histNMC
1523	#include "ini_histhfNMC.h"
1524	#include "ini_histdayNMC.h"
1525	#include "ini_histmthNMC.h"
1526	#endif
1527
1528	#include "ini_histday_seri.h"
1529
1530	#include "ini_paramLMDZ_phy.h"
1531
1532	#endif
1533
1534	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1535	ecrit_hf2mth = ecrit_mth/ecrit_hf
1536
1537	ecrit_hf = ecrit_hf * un_jour
1538	cIM
1539	IF(ecrit_day.LE.1.) THEN
1540	ecrit_day = ecrit_day * un_jour !en secondes
1541	ENDIF
1542	cIM
1543	ecrit_mth = ecrit_mth * un_jour
1544	ecrit_ins = ecrit_ins * un_jour
1545	ecrit_reg = ecrit_reg * un_jour
1546	ecrit_tra = ecrit_tra * un_jour
1547	ecrit_ISCCP = ecrit_ISCCP * un_jour
1548	ecrit_LES = ecrit_LES * un_jour
1549	c
1550	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1551	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1552	. ecrit_hf2mth
1553	cIM 030306 END
1554
1555
1556	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1557	date0 = jD_ref
1558	WRITE(,) 'physiq date0 : ',date0
1559	c
1560	c
1561	c
1562	c Prescrire l'ozone dans l'atmosphere
1563	c
1564	c
1565	cc DO i = 1, klon
1566	cc DO k = 1, klev
1567	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1568	cc ENDDO
1569	cc ENDDO
1570	c
1571	IF (config_inca /= 'none') THEN
1572	#ifdef INCA
1573	CALL VTe(VTphysiq)
1574	CALL VTb(VTinca)
1575	! iii = MOD(NINT(xjour),360)
1576	! calday = REAL(iii) + jH_cur
1577	calday = REAL(days_elapsed) + jH_cur
1578	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1579
1580	CALL chemini(
1581	$ rg,
1582	$ ra,
1583	$ airephy,
1584	$ rlat,
1585	$ rlon,
1586	$ presnivs,
1587	$ calday,
1588	$ klon,
1589	$ nqtot,
1590	$ pdtphys,
1591	$ annee_ref,
1592	$ day_ref,
1593	$ itau_phy)
1594
1595	CALL VTe(VTinca)
1596	CALL VTb(VTphysiq)
1597	#endif
1598	END IF
1599	c
1600	c
1601	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1602	! Nouvelle initialisation pour le rayonnement RRTM
1603	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1604
1605	call iniradia(klon,klev,paprs(1,1:klev+1))
1606
1607	C$omp single
1608	if (read_climoz >= 1) then
1609	call open_climoz(ncid_climoz, press_climoz)
1610	END IF
1611	C$omp end single
1612	ENDIF
1613	!
1614	! ************** Fin de IF ( debut ) *************
1615	!
1616	!
1617	! Incrementer le compteur de la physique
1618	!
1619	itap = itap + 1
1620	!
1621	! Update fraction of the sub-surfaces (pctsrf) and
1622	! initialize, where a new fraction has appeared, all variables depending
1623	! on the surface fraction.
1624	!
1625	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1626	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1627
1628	! Tendances bidons pour les processus qui n'affectent pas certaines
1629	! variables.
1630	du0(:,:)=0.
1631	dv0(:,:)=0.
1632	dq0(:,:)=0.
1633	dql0(:,:)=0.
1634	c
1635	c Mettre a zero des variables de sortie (pour securite)
1636	c
1637	DO i = 1, klon
1638	d_ps(i) = 0.0
1639	ENDDO
1640	DO k = 1, klev
1641	DO i = 1, klon
1642	d_t(i,k) = 0.0
1643	d_u(i,k) = 0.0
1644	d_v(i,k) = 0.0
1645	ENDDO
1646	ENDDO
1647	DO iq = 1, nqtot
1648	DO k = 1, klev
1649	DO i = 1, klon
1650	d_qx(i,k,iq) = 0.0
1651	ENDDO
1652	ENDDO
1653	ENDDO
1654	da(:,:)=0.
1655	mp(:,:)=0.
1656	phi(:,:,:)=0.
1657	c
1658	c Ne pas affecter les valeurs entrees de u, v, h, et q
1659	c
1660	DO k = 1, klev
1661	DO i = 1, klon
1662	t_seri(i,k) = t(i,k)
1663	u_seri(i,k) = u(i,k)
1664	v_seri(i,k) = v(i,k)
1665	q_seri(i,k) = qx(i,k,ivap)
1666	ql_seri(i,k) = qx(i,k,iliq)
1667	qs_seri(i,k) = 0.
1668	ENDDO
1669	ENDDO
1670	IF (nqtot.GE.3) THEN
1671	DO iq = 3, nqtot
1672	DO k = 1, klev
1673	DO i = 1, klon
1674	tr_seri(i,k,iq-2) = qx(i,k,iq)
1675	ENDDO
1676	ENDDO
1677	ENDDO
1678	ELSE
1679	DO k = 1, klev
1680	DO i = 1, klon
1681	tr_seri(i,k,1) = 0.0
1682	ENDDO
1683	ENDDO
1684	ENDIF
1685	C
1686	DO i = 1, klon
1687	ztsol(i) = 0.
1688	ENDDO
1689	DO nsrf = 1, nbsrf
1690	DO i = 1, klon
1691	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1692	ENDDO
1693	ENDDO
1694	cIM
1695	IF (ip_ebil_phy.ge.1) THEN
1696	ztit='after dynamic'
1697	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1698	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1699	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1700	C Comme les tendances de la physique sont ajoute dans la dynamique,
1701	C on devrait avoir que la variation d'entalpie par la dynamique
1702	C est egale a la variation de la physique au pas de temps precedent.
1703	C Donc la somme de ces 2 variations devrait etre nulle.
1704	call diagphy(airephy,ztit,ip_ebil_phy
1705	e , zero_v, zero_v, zero_v, zero_v, zero_v
1706	e , zero_v, zero_v, zero_v, ztsol
1707	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1708	s , fs_bound, fq_bound )
1709	END IF
1710
1711	c Diagnostiquer la tendance dynamique
1712	c
1713	IF (ancien_ok) THEN
1714	DO k = 1, klev
1715	DO i = 1, klon
1716	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1717	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1718	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1719	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1720	ENDDO
1721	ENDDO
1722	ELSE
1723	DO k = 1, klev
1724	DO i = 1, klon
1725	d_u_dyn(i,k) = 0.0
1726	d_v_dyn(i,k) = 0.0
1727	d_t_dyn(i,k) = 0.0
1728	d_q_dyn(i,k) = 0.0
1729	ENDDO
1730	ENDDO
1731	ancien_ok = .TRUE.
1732	ENDIF
1733	c
1734	c Ajouter le geopotentiel du sol:
1735	c
1736	DO k = 1, klev
1737	DO i = 1, klon
1738	zphi(i,k) = pphi(i,k) + pphis(i)
1739	ENDDO
1740	ENDDO
1741	c
1742	c Verifier les temperatures
1743	c
1744	cIM BEG
1745	IF (check) THEN
1746	amn=MIN(ftsol(1,is_ter),1000.)
1747	amx=MAX(ftsol(1,is_ter),-1000.)
1748	DO i=2, klon
1749	amn=MIN(ftsol(i,is_ter),amn)
1750	amx=MAX(ftsol(i,is_ter),amx)
1751	ENDDO
1752	c
1753	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1754	ENDIF !(check) THEN
1755	cIM END
1756	c
1757	CALL hgardfou(t_seri,ftsol,'debutphy')
1758	c
1759	cIM BEG
1760	IF (check) THEN
1761	amn=MIN(ftsol(1,is_ter),1000.)
1762	amx=MAX(ftsol(1,is_ter),-1000.)
1763	DO i=2, klon
1764	amn=MIN(ftsol(i,is_ter),amn)
1765	amx=MAX(ftsol(i,is_ter),amx)
1766	ENDDO
1767	c
1768	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1769	ENDIF !(check) THEN
1770	cIM END
1771	c
1772	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1773	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1774	c
1775	if (read_climoz >= 1) then
1776	C Ozone from a file
1777	! Update required ozone index:
1778	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1779	$ / ioget_year_len(year_cur) * 360.) + 1
1780	if (ro3i == 361) ro3i = 360
1781	C (This should never occur, except perhaps because of roundup
1782	C error. See documentation.)
1783	if (ro3i /= co3i) then
1784	C Update ozone field:
1785	if (read_climoz == 1) then
1786	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1787	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1788	else
1789	C read_climoz == 2
1790	call regr_pr_av(ncid_climoz,
1791	$ (/"tro3 ", "tro3_daylight"/),
1792	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1793	$ v3=wo)
1794	end if
1795	! Convert from mole fraction of ozone to column density of ozone in a
1796	! cell, in kDU:
1797	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1798	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1799	C (By regridding ozone values for LMDZ only once every 360th of
1800	C year, we have already neglected the variation of pressure in one
1801	C 360th of year. So do not recompute "wo" at each time step even if
1802	C "zmasse" changes a little.)
1803	co3i = ro3i
1804	end if
1805	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1806	C Once per day, update ozone from Royer:
1807	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1808	ENDIF
1809	c
1810	c Re-evaporer l'eau liquide nuageuse
1811	c
1812	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1813	DO i = 1, klon
1814	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1815	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1816	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1817	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1818	zb = MAX(0.0,ql_seri(i,k))
1819	za = - MAX(0.0,ql_seri(i,k))
1820	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1821	t_seri(i,k) = t_seri(i,k) + za
1822	q_seri(i,k) = q_seri(i,k) + zb
1823	ql_seri(i,k) = 0.0
1824	d_t_eva(i,k) = za
1825	d_q_eva(i,k) = zb
1826	ENDDO
1827	ENDDO
1828	cIM
1829	IF (ip_ebil_phy.ge.2) THEN
1830	ztit='after reevap'
1831	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1832	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1833	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1834	call diagphy(airephy,ztit,ip_ebil_phy
1835	e , zero_v, zero_v, zero_v, zero_v, zero_v
1836	e , zero_v, zero_v, zero_v, ztsol
1837	e , d_h_vcol, d_qt, d_ec
1838	s , fs_bound, fq_bound )
1839	C
1840	END IF
1841
1842	c
1843	c=========================================================================
1844	! Calculs de l'orbite.
1845	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1846	! doit donc etre placé avant radlwsw et pbl_surface
1847
1848	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1849	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1850	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1851	!
1852	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1853	! solarlong0
1854	if (solarlong0<-999.) then
1855	if (new_orbit) then
1856	! calcul selon la routine utilisee pour les planetes
1857	call solarlong(day_since_equinox, zlongi, dist)
1858	else
1859	! calcul selon la routine utilisee pour l'AR4
1860	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1861	endif
1862	else
1863	zlongi=solarlong0 ! longitude solaire vraie
1864	dist=1. ! distance au soleil / moyenne
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	I rlon,rlat,qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3380
3381
3382	ENDIF
3383
3384	c
3385	c Calculer le transport de l'eau et de l'energie (diagnostique)
3386	c
3387	CALL transp (paprs,zxtsol,
3388	e t_seri, q_seri, u_seri, v_seri, zphi,
3389	s ve, vq, ue, uq)
3390	c
3391	cIM global posePB BEG
3392	IF(1.EQ.0) THEN
3393	c
3394	CALL transp_lay (paprs,zxtsol,
3395	e t_seri, q_seri, u_seri, v_seri, zphi,
3396	s ve_lay, vq_lay, ue_lay, uq_lay)
3397	c
3398	ENDIF !(1.EQ.0) THEN
3399	cIM global posePB END
3400	c Accumuler les variables a stocker dans les fichiers histoire:
3401	c
3402	c+jld ec_conser
3403	DO k = 1, klev
3404	DO i = 1, klon
3405	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3406	d_t_ec(i,k)=0.5/ZRCPD
3407	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3408	ENDDO
3409	ENDDO
3410
3411	DO k = 1, klev
3412	DO i = 1, klon
3413	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3414	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3415	END DO
3416	END DO
3417	c-jld ec_conser
3418	cIM
3419	IF (ip_ebil_phy.ge.1) THEN
3420	ztit='after physic'
3421	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3422	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3423	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3424	C Comme les tendances de la physique sont ajoute dans la dynamique,
3425	C on devrait avoir que la variation d'entalpie par la dynamique
3426	C est egale a la variation de la physique au pas de temps precedent.
3427	C Donc la somme de ces 2 variations devrait etre nulle.
3428
3429	call diagphy(airephy,ztit,ip_ebil_phy
3430	e , topsw, toplw, solsw, sollw, sens
3431	e , evap, rain_fall, snow_fall, ztsol
3432	e , d_h_vcol, d_qt, d_ec
3433	s , fs_bound, fq_bound )
3434	C
3435	d_h_vcol_phy=d_h_vcol
3436	C
3437	END IF
3438	C
3439	c=======================================================================
3440	c SORTIES
3441	c=======================================================================
3442
3443	cIM Interpolation sur les niveaux de pression du NMC
3444	c -------------------------------------------------
3445	c
3446	#include "calcul_STDlev.h"
3447	twriteSTD(:,:,1)=tsumSTD(:,:,1)
3448	qwriteSTD(:,:,1)=qsumSTD(:,:,1)
3449	rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
3450	phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
3451	uwriteSTD(:,:,1)=usumSTD(:,:,1)
3452	vwriteSTD(:,:,1)=vsumSTD(:,:,1)
3453	wwriteSTD(:,:,1)=wsumSTD(:,:,1)
3454
3455	twriteSTD(:,:,2)=tsumSTD(:,:,2)
3456	qwriteSTD(:,:,2)=qsumSTD(:,:,2)
3457	rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
3458	phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
3459	uwriteSTD(:,:,2)=usumSTD(:,:,2)
3460	vwriteSTD(:,:,2)=vsumSTD(:,:,2)
3461	wwriteSTD(:,:,2)=wsumSTD(:,:,2)
3462
3463	twriteSTD(:,:,3)=tlevSTD(:,:)
3464	qwriteSTD(:,:,3)=qlevSTD(:,:)
3465	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3466	phiwriteSTD(:,:,3)=philevSTD(:,:)
3467	uwriteSTD(:,:,3)=ulevSTD(:,:)
3468	vwriteSTD(:,:,3)=vlevSTD(:,:)
3469	wwriteSTD(:,:,3)=wlevSTD(:,:)
3470
3471	twriteSTD(:,:,4)=tlevSTD(:,:)
3472	qwriteSTD(:,:,4)=qlevSTD(:,:)
3473	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3474	phiwriteSTD(:,:,4)=philevSTD(:,:)
3475	uwriteSTD(:,:,4)=ulevSTD(:,:)
3476	vwriteSTD(:,:,4)=vlevSTD(:,:)
3477	wwriteSTD(:,:,4)=wlevSTD(:,:)
3478	c
3479	cIM initialisation 5eme fichier de sortie
3480	cIM ajoute 5eme niveau 170310 BEG
3481	twriteSTD(:,:,5)=tlevSTD(:,:)
3482	qwriteSTD(:,:,5)=qlevSTD(:,:)
3483	rhwriteSTD(:,:,5)=rhlevSTD(:,:)
3484	phiwriteSTD(:,:,5)=philevSTD(:,:)
3485	uwriteSTD(:,:,5)=ulevSTD(:,:)
3486	vwriteSTD(:,:,5)=vlevSTD(:,:)
3487	wwriteSTD(:,:,5)=wlevSTD(:,:)
3488	cIM for NMC files
3489	DO n=1, nlevSTD3
3490	DO k=1, nlevSTD
3491	if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
3492	twriteSTD3(:,n)=tlevSTD(:,k)
3493	qwriteSTD3(:,n)=qlevSTD(:,k)
3494	rhwriteSTD3(:,n)=rhlevSTD(:,k)
3495	phiwriteSTD3(:,n)=philevSTD(:,k)
3496	uwriteSTD3(:,n)=ulevSTD(:,k)
3497	vwriteSTD3(:,n)=vlevSTD(:,k)
3498	wwriteSTD3(:,n)=wlevSTD(:,k)
3499	endif !rlevSTD3(n).EQ.rlevSTD(k)
3500	ENDDO
3501	ENDDO
3502	c
3503	DO n=1, nlevSTD8
3504	DO k=1, nlevSTD
3505	if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
3506	tnondefSTD8(:,n)=tnondef(:,k,2)
3507	twriteSTD8(:,n)=tsumSTD(:,k,2)
3508	qwriteSTD8(:,n)=qsumSTD(:,k,2)
3509	rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
3510	phiwriteSTD8(:,n)=phisumSTD(:,k,2)
3511	uwriteSTD8(:,n)=usumSTD(:,k,2)
3512	vwriteSTD8(:,n)=vsumSTD(:,k,2)
3513	wwriteSTD8(:,n)=wsumSTD(:,k,2)
3514	endif !rlevSTD8(n).EQ.rlevSTD(k)
3515	ENDDO
3516	ENDDO
3517	c
3518	c slp sea level pressure
3519	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3520	c
3521	ccc prw = eau precipitable
3522	DO i = 1, klon
3523	prw(i) = 0.
3524	DO k = 1, klev
3525	prw(i) = prw(i) +
3526	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3527	ENDDO
3528	ENDDO
3529	c
3530	cIM initialisation + calculs divers diag AMIP2
3531	c
3532	#include "calcul_divers.h"
3533	c
3534	IF (config_inca /= 'none') THEN
3535	#ifdef INCA
3536	CALL VTe(VTphysiq)
3537	CALL VTb(VTinca)
3538
3539	CALL chemhook_end (
3540	$ dtime,
3541	$ pplay,
3542	$ t_seri,
3543	$ tr_seri,
3544	$ nbtr,
3545	$ paprs,
3546	$ q_seri,
3547	$ airephy,
3548	$ pphi,
3549	$ pphis,
3550	$ zx_rh)
3551
3552	CALL VTe(VTinca)
3553	CALL VTb(VTphysiq)
3554	#endif
3555	END IF
3556
3557	c=============================================================
3558	c
3559	c Convertir les incrementations en tendances
3560	c
3561	if (mydebug) then
3562	call writefield_phy('u_seri',u_seri,llm)
3563	call writefield_phy('v_seri',v_seri,llm)
3564	call writefield_phy('t_seri',t_seri,llm)
3565	call writefield_phy('q_seri',q_seri,llm)
3566	endif
3567
3568	DO k = 1, klev
3569	DO i = 1, klon
3570	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3571	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3572	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3573	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3574	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3575	ENDDO
3576	ENDDO
3577	c
3578	IF (nqtot.GE.3) THEN
3579	DO iq = 3, nqtot
3580	DO k = 1, klev
3581	DO i = 1, klon
3582	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3583	ENDDO
3584	ENDDO
3585	ENDDO
3586	ENDIF
3587	c
3588	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3589	cIM global posePB#include "write_bilKP_ins.h"
3590	cIM global posePB#include "write_bilKP_ave.h"
3591	c
3592
3593	c Sauvegarder les valeurs de t et q a la fin de la physique:
3594	c
3595	DO k = 1, klev
3596	DO i = 1, klon
3597	u_ancien(i,k) = u_seri(i,k)
3598	v_ancien(i,k) = v_seri(i,k)
3599	t_ancien(i,k) = t_seri(i,k)
3600	q_ancien(i,k) = q_seri(i,k)
3601	ENDDO
3602	ENDDO
3603	c
3604	!==========================================================================
3605	! Sorties des tendances pour un point particulier
3606	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3607	! pour le debug
3608	! La valeur de igout est attribuee plus haut dans le programme
3609	!==========================================================================
3610
3611	if (prt_level.ge.1) then
3612	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3613	write(lunout,*)
3614	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3615	write(lunout,*)
3616	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3617	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3618	s pctsrf(igout,is_sic)
3619	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3620	do k=1,klev
3621	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3622	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3623	s d_t_eva(igout,k)
3624	enddo
3625	write(lunout,*) 'cool,heat'
3626	do k=1,klev
3627	write(lunout,*) cool(igout,k),heat(igout,k)
3628	enddo
3629
3630	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3631	do k=1,klev
3632	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3633	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3634	enddo
3635
3636	write(lunout,*) 'd_ps ',d_ps(igout)
3637	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3638	do k=1,klev
3639	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3640	s d_qx(igout,k,1),d_qx(igout,k,2)
3641	enddo
3642	endif
3643
3644	!==========================================================================
3645
3646	c============================================================
3647	c Calcul de la temperature potentielle
3648	c============================================================
3649	DO k = 1, klev
3650	DO i = 1, klon
3651	cJYG/IM theta en debut du pas de temps
3652	cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3653	cJYG/IM theta en fin de pas de temps de physique
3654	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
3655	ENDDO
3656	ENDDO
3657	c
3658
3659	c 22.03.04 BEG
3660	c=============================================================
3661	c Ecriture des sorties
3662	c=============================================================
3663	#ifdef CPP_IOIPSL
3664
3665	c Recupere des varibles calcule dans differents modules
3666	c pour ecriture dans histxxx.nc
3667
3668	! Get some variables from module fonte_neige_mod
3669	CALL fonte_neige_get_vars(pctsrf,
3670	. zxfqcalving, zxfqfonte, zxffonte)
3671
3672
3673	#include "phys_output_write.h"
3674
3675	#ifdef histISCCP
3676	#include "write_histISCCP.h"
3677	#endif
3678
3679	#ifdef histNMC
3680	#include "write_histhfNMC.h"
3681	#include "write_histdayNMC.h"
3682	#include "write_histmthNMC.h"
3683	#endif
3684
3685	#include "write_histday_seri.h"
3686
3687	#include "write_paramLMDZ_phy.h"
3688
3689	#endif
3690
3691	c 22.03.04 END
3692	c
3693	c====================================================================
3694	c Si c'est la fin, il faut conserver l'etat de redemarrage
3695	c====================================================================
3696	c
3697
3698
3699	IF (lafin) THEN
3700	itau_phy = itau_phy + itap
3701	CALL phyredem ("restartphy.nc")
3702	! open(97,form="unformatted",file="finbin")
3703	! write(97) u_seri,v_seri,t_seri,q_seri
3704	! close(97)
3705	C$OMP MASTER
3706	if (read_climoz >= 1) then
3707	if (is_mpi_root) then
3708	call nf95_close(ncid_climoz)
3709	end if
3710	deallocate(press_climoz) ! pointer
3711	end if
3712	C$OMP END MASTER
3713	ENDIF
3714
3715	! first=.false.
3716
3717	RETURN
3718	END
3719	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3720	IMPLICIT none
3721	c
3722	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3723	c la conservation de l'eau
3724	c
3725	#include "YOMCST.h"
3726	INTEGER klon,klev
3727	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3728	REAL aire(klon)
3729	REAL qtotal, zx, qcheck
3730	INTEGER i, k
3731	c
3732	zx = 0.0
3733	DO i = 1, klon
3734	zx = zx + aire(i)
3735	ENDDO
3736	qtotal = 0.0
3737	DO k = 1, klev
3738	DO i = 1, klon
3739	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3740	. *(paprs(i,k)-paprs(i,k+1))/RG
3741	ENDDO
3742	ENDDO
3743	c
3744	qcheck = qtotal/zx
3745	c
3746	RETURN
3747	END
3748	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3749	IMPLICIT none
3750	c
3751	c Tranformer une variable de la grille physique a
3752	c la grille d'ecriture
3753	c
3754	INTEGER nfield,nlon,iim,jjmp1, jjm
3755	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3756	c
3757	INTEGER i, n, ig
3758	c
3759	jjm = jjmp1 - 1
3760	DO n = 1, nfield
3761	DO i=1,iim
3762	ecrit(i,n) = fi(1,n)
3763	ecrit(i+jjm*iim,n) = fi(nlon,n)
3764	ENDDO
3765	DO ig = 1, nlon - 2
3766	ecrit(iim+ig,n) = fi(1+ig,n)
3767	ENDDO
3768	ENDDO
3769	RETURN
3770	END

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

Download in other formats:

Original Format