Context Navigation

physiq.F @ 1507

Last change on this file since 1507 was 1507, checked in by idelkadi, 13 years ago

Introduction de nouvelles facon de calucler les ratqs pour firstilp quand iflag_cldcon>=7 : ratqs devient pronos

tic avec une relaxation vers un ratqs issu de la convection (ratqsc) quand celle-ci est active et une relaxation ve
rs ratqss dans le cas contraire.
Les constantes de temps sont a 3h en dur dans le code (en attendant mieux)
On met un coeff iflag_cldcon/100 devant ratqsc
On plafone le ratqs a 0.5 (la aussi en dur dans le code en attendant de stabiliser une facon de faire).

Introduction de nouveaux diagnostics concernant la separation entre la contribution des thermiques et le reste d

ans fisrtilp (dq/tlscth, dq/tlscst,plulth, plulst)
Correction de la sortie de la hauteur des thermiques zmax_th

Introduction d'une option iflag_wake=2 pour prendre en compte les precipitations issues de fisrtipl dans le forc

age des wake.

Correction d'un bug sur l'identification des points affectes par les thermiques (ptconvth(:,k)=fm_therm(:,k+1)>0

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

File size: 128.5 KB

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

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

Download in other formats:

Original Format