Context Navigation

physiq.F @ 1496

Last change on this file since 1496 was 1496, checked in by Laurent Fairhead, 13 years ago

Modified closure for the new physics package, new values for the iflag_coupl parameter
FH

Fermeture modifiée pour la nouvelle physique, nouvelles valeurs définies pour
le paramètre iflag_coupl
FH

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

File size: 126.3 KB

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

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

Download in other formats:

Original Format