Context Navigation

physiq.F90 @ 151

Last change on this file since 151 was 151, checked in by lfita, 10 years ago
Increasing value of 'largest'
File size: 182.1 KB

Line
1	! $Id: physiq.F 1795 2013-07-18 08:20:28Z emillour $
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	!L. Fita, LMD. November 2013
26	! USE pbl_surface_mod, ONLY : pbl_surface
27	! USE pbl_surface_mod, ONLY : pbl_surface, pbl_surface_init
28	USE pbl_surface_mod
29	USE change_srf_frac_mod
30	USE surface_data, ONLY : type_ocean, ok_veget
31	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
32	USE phys_state_var_mod ! Variables sauvegardees de la physique
33	USE phys_output_var_mod ! Variables pour les ecritures des sorties
34	USE fonte_neige_mod, ONLY : fonte_neige_get_vars, fonte_neige_init
35	USE phys_output_mod
36	USE phys_output_ctrlout_mod
37	USE iophy
38	use open_climoz_m, only: open_climoz ! ozone climatology from a file
39	use regr_pr_av_m, only: regr_pr_av
40	use netcdf95, only: nf95_close
41	!IM for NMC files
42	!c use netcdf, only: nf90_fill_real
43	use netcdf
44	use mod_phys_lmdz_mpi_data, only: is_mpi_root
45	USE aero_mod
46	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
47	use conf_phys_m, only: conf_phys
48	use radlwsw_m, only: radlwsw
49	USE control_mod
50	#ifdef REPROBUS
51	USE CHEM_REP, ONLY : Init_chem_rep_xjour
52	#endif
53	USE indice_sol_mod
54
55	!IM stations CFMIP
56	USE CFMIP_point_locations
57
58	!L. Fita. July 2013
59	USE orografi
60	USE orografi_strato_mod
61	USE diagphy_mod
62	! USE orbite_mod, ONLY : angle
63	USE phyaqua_mod, ONLY: zenang_an
64	USE nuage_mod, ONLY : diagcld1
65
66	! L. Fita, LMD. November 2013.
67	!! USE wrf_lmdz_mod
68	!! USE module_lmdz_variables, ONLY: neige_initialize, limit_initialize
69	!! USE physiq_limit_variables_mod
70	!! USE physiq_limit_variables_mod
71	USE lmdz_wrf_variables_mod
72	! L. Fita, LMD. January 2014
73	! Defining with SAVE statement all that variables that should appear in the output
74	USE output_lmdz_NOmodule
75
76	IMPLICIT none
77	!>======================================================================
78	!!
79	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
80	!!
81	!! Objet: Moniteur general de la physique du modele
82	!!AA Modifications quant aux traceurs :
83	!!AA - uniformisation des parametrisations ds phytrac
84	!!AA - stockage des moyennes des champs necessaires
85	!!AA en mode traceur off-line
86	!!======================================================================
87	!! CLEFS CPP POUR LES IO
88	!! =====================
89	#define histNMC
90	!c#define histISCCP
91	!!======================================================================
92	!! modif ( P. Le Van , 12/10/98 )
93	!!
94	!! Arguments:
95	!!
96	!! nlon----input-I-nombre de points horizontaux
97	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
98	!! debut---input-L-variable logique indiquant le premier passage
99	!! lafin---input-L-variable logique indiquant le dernier passage
100	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
101	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
102	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
103	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
104	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
105	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
106	!! pphis---input-R-geopotentiel du sol
107	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
108	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
109	!! v-------input-R-vitesse Y (de S a N) en m/s
110	!! t-------input-R-temperature (K)
111	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
112	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
113	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
114	!! flxmass_w -input-R- flux de masse verticale
115	!! d_u-----output-R-tendance physique de "u" (m/s/s)
116	!! d_v-----output-R-tendance physique de "v" (m/s/s)
117	!! d_t-----output-R-tendance physique de "t" (K/s)
118	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
119	!! d_ps----output-R-tendance physique de la pression au sol
120	!!IM
121	!! PVteta--output-R-vorticite potentielle a des thetas constantes
122	!!======================================================================
123	#include "dimensions.h"
124	integer jjmp1
125	parameter (jjmp1=jjm+1-1/jjm)
126	integer iip1
127	parameter (iip1=iim+1)
128
129	#include "regdim.h"
130	#include "dimsoil.h"
131	#include "clesphys.h"
132	#include "temps.h"
133	#include "iniprint.h"
134	#include "thermcell.h"
135	!c======================================================================
136	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
137	PARAMETER (ok_cvl=.TRUE.)
138	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
139	PARAMETER (ok_gust=.FALSE.)
140	integer iflag_radia ! active ou non le rayonnement (MPL)
141	save iflag_radia
142	!$OMP THREADPRIVATE(iflag_radia)
143	!c======================================================================
144	LOGICAL check ! Verifier la conservation du modele en eau
145	PARAMETER (check=.FALSE.)
146	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
147	PARAMETER (ok_stratus=.FALSE.)
148	!c======================================================================
149	REAL amn, amx
150	INTEGER igout
151	!c======================================================================
152	!c Clef controlant l'activation du cycle diurne:
153	!ccc LOGICAL cycle_diurne
154	!ccc PARAMETER (cycle_diurne=.FALSE.)
155	!c======================================================================
156	!c Modele thermique du sol, a activer pour le cycle diurne:
157	!ccc LOGICAL soil_model
158	!ccc PARAMETER (soil_model=.FALSE.)
159	!c======================================================================
160	!c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
161	!c le calcul du rayonnement est celle apres la precipitation des nuages.
162	!c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
163	!c la condensation et la precipitation. Cette cle augmente les impacts
164	!c radiatifs des nuages.
165	!ccc LOGICAL new_oliq
166	!ccc PARAMETER (new_oliq=.FALSE.)
167	!c======================================================================
168	!c Clefs controlant deux parametrisations de l'orographie:
169	!cc LOGICAL ok_orodr
170	!ccc PARAMETER (ok_orodr=.FALSE.)
171	!ccc LOGICAL ok_orolf
172	!ccc PARAMETER (ok_orolf=.FALSE.)
173	!c======================================================================
174	LOGICAL ok_journe ! sortir le fichier journalier
175	save ok_journe
176	!$OMP THREADPRIVATE(ok_journe)
177	!c
178	LOGICAL ok_mensuel ! sortir le fichier mensuel
179	save ok_mensuel
180	!$OMP THREADPRIVATE(ok_mensuel)
181	!c
182	LOGICAL ok_instan ! sortir le fichier instantane
183	save ok_instan
184	!$OMP THREADPRIVATE(ok_instan)
185	!c
186	LOGICAL ok_LES ! sortir le fichier LES
187	save ok_LES
188	!$OMP THREADPRIVATE(ok_LES)
189	!c
190	LOGICAL callstats ! sortir le fichier stats
191	save callstats
192	!$OMP THREADPRIVATE(callstats)
193	!c
194	LOGICAL ok_region ! sortir le fichier regional
195	PARAMETER (ok_region=.FALSE.)
196	!c======================================================================
197	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
198	! real weak_inversion(klon),dthmin(klon)
199	real seuil_inversion
200	save seuil_inversion
201	!$OMP THREADPRIVATE(seuil_inversion)
202	integer iflag_ratqs
203	save iflag_ratqs
204	!$OMP THREADPRIVATE(iflag_ratqs)
205	real facteur
206
207	REAL zz,znum,zden
208	REAL wmax_th(klon)
209	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
210	! REAL zmax_th(klon)
211	REAL tau_overturning_th(klon)
212
213	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
214	! integer lmax_th(klon)
215	integer limbas(klon)
216	real ratqscth(klon,klev)
217	real ratqsdiff(klon,klev)
218	real zqsatth(klon,klev)
219
220	!c======================================================================
221	!c
222	INTEGER ivap ! indice de traceurs pour vapeur d'eau
223	PARAMETER (ivap=1)
224	INTEGER iliq ! indice de traceurs pour eau liquide
225	PARAMETER (iliq=2)
226
227	!c
228	!c
229	!c Variables argument:
230	!c
231	INTEGER nlon
232	INTEGER nlev
233	REAL, intent(in):: jD_cur, jH_cur
234
235	REAL pdtphys
236	LOGICAL debut, lafin
237	REAL paprs(klon,klev+1)
238	REAL pplay(klon,klev)
239	REAL pphi(klon,klev)
240	REAL pphis(klon)
241	REAL presnivs(klev)
242	REAL znivsig(klev)
243	real pir
244
245	REAL u(klon,klev)
246	REAL v(klon,klev)
247	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
248	! REAL t(klon,klev),theta(klon,klev),thetal(klon,klev)
249	REAL t(klon,klev),thetal(klon,klev)
250
251	!c thetal: ligne suivante a decommenter si vous avez les fichiers MPL 20130625
252	!c fth_fonctions.F90 et parkind1.F90
253	!c sinon thetal=theta
254	!c REAL fth_thetae,fth_thetav,fth_thetal
255	REAL qx(klon,klev,nqtot)
256	REAL flxmass_w(klon,klev)
257	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
258	! REAL omega(klon,klev) ! vitesse verticale en Pa/s
259	REAL d_u(klon,klev)
260	REAL d_v(klon,klev)
261	REAL d_t(klon,klev)
262	REAL d_qx(klon,klev,nqtot)
263	REAL d_ps(klon)
264	! Variables pour le transport convectif
265	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
266	! Variables pour le lessivage convectif
267	! RomP >>>
268	real phi2(klon,klev,klev)
269	real d1a(klon,klev),dam(klon,klev)
270	real ev(klon,klev),ep(klon,klev)
271	real clw(klon,klev),elij(klon,klev,klev)
272	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
273	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
274	! real wdtrainA(klon,klev),wdtrainM(klon,klev)
275
276	! RomP <<<
277	!IM definition dynamique o_trac dans phys_output_open
278	! type(ctrl_out) :: o_trac(nqtot)
279	!c
280	!IM Amip2 PV a theta constante
281	!c
282	INTEGER nbteta
283	PARAMETER(nbteta=3)
284	CHARACTER*3 ctetaSTD(nbteta)
285	DATA ctetaSTD/'350','380','405'/
286	SAVE ctetaSTD
287	!$OMP THREADPRIVATE(ctetaSTD)
288	REAL rtetaSTD(nbteta)
289	DATA rtetaSTD/350., 380., 405./
290	SAVE rtetaSTD
291	!$OMP THREADPRIVATE(rtetaSTD)
292	!c
293	REAL PVteta(klon,nbteta)
294	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
295	!c
296	!MI Amip2 PV a theta constante
297
298	!cym INTEGER klevp1, klevm1
299	!cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
300	!cym#include "raddim.h"
301	!c
302	!c
303	!IM Amip2
304	!c variables a une pression donnee
305	!c
306	#include "declare_STDlev.h"
307	!c
308	CHARACTER*4 bb2
309	CHARACTER*2 bb3
310	!c
311	#include "radopt.h"
312	!c
313	!c
314	!c prw: precipitable water
315	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
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	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
330	! REAL flwp(klon), fiwp(klon)
331	! REAL flwc(klon,klev), fiwc(klon,klev)
332	REAL flwp_c(klon), fiwp_c(klon)
333	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
334	REAL flwp_s(klon), fiwp_s(klon)
335	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
336
337	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
338	! REAL evap_pot(klon,nbsrf)
339
340	!IM ISCCP simulator v3.4
341	!c dans clesphys.h top_height, overlap
342	!cv3.4
343	INTEGER debug, debugcol
344	!cym INTEGER npoints
345	!cym PARAMETER(npoints=klon)
346	!c
347	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
348	INTEGER nregISCtot
349	PARAMETER(nregISCtot=1)
350	!c
351	!c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
352	!c y compris pour 1 point
353	!c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
354	!c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
355	INTEGER imin_debut, nbpti
356	INTEGER jmin_debut, nbptj
357	!IM parametres ISCCP BEG
358	INTEGER nbapp_isccp
359	! INTEGER nbapp_isccp,isccppas
360	! PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
361	! !i.e. toutes les 3 heures
362	INTEGER n
363	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
364	DATA ifreq_isccp/3/
365	SAVE ifreq_isccp
366	!$OMP THREADPRIVATE(ifreq_isccp)
367	CHARACTER*5 typinout(napisccp)
368	DATA typinout/'i3od'/
369	SAVE typinout
370	!$OMP THREADPRIVATE(typinout)
371	!IM verif boxptop BEG
372	CHARACTER*1 verticaxe(napisccp)
373	DATA verticaxe/'1'/
374	SAVE verticaxe
375	!$OMP THREADPRIVATE(verticaxe)
376	!IM verif boxptop END
377	INTEGER nvlev(napisccp)
378	!c INTEGER nvlev
379	REAL t1, aa
380	REAL seed_re(klon,napisccp)
381	!cym !!!! A voir plus tard
382	!cym INTEGER iphy(iim,jjmp1)
383	!IM parametres ISCCP END
384	!c
385	!c ncol = nb. de sous-colonnes pour chaque maille du GCM
386	!c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
387	!c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
388	INTEGER,SAVE :: ncol(napisccp)
389	!$OMP THREADPRIVATE(ncol)
390	INTEGER ncolmx, seed(klon,napisccp)
391	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
392	!c PARAMETER(ncolmx=1500)
393	PARAMETER(ncolmx=300)
394	!c
395	!IM verif boxptop BEG
396	REAL vertlev(ncolmx,napisccp)
397	!IM verif boxptop END
398	!c
399	REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
400	INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
401	!$OMP THREADPRIVATE(tautab,invtau)
402	REAL emsfc_lw
403	PARAMETER(emsfc_lw=0.99)
404	!c REAL ran0 ! type for random number fuction
405	!c
406	REAL cldtot(klon,klev)
407	!c variables de haut en bas pour le simulateur ISCCP
408	REAL dtau_s(klon,klev) !tau nuages startiformes
409	REAL dtau_c(klon,klev) !tau nuages convectifs
410	REAL dem_s(klon,klev) !emissivite nuages startiformes
411	REAL dem_c(klon,klev) !emissivite nuages convectifs
412	!c
413	!c variables de haut en bas pour le simulateur ISCCP
414	REAL pfull(klon,klev)
415	REAL phalf(klon,klev+1)
416	REAL qv(klon,klev)
417	REAL cc(klon,klev)
418	REAL conv(klon,klev)
419	REAL dtau_sH2B(klon,klev)
420	REAL dtau_cH2B(klon,klev)
421	REAL at(klon,klev)
422	REAL dem_sH2B(klon,klev)
423	REAL dem_cH2B(klon,klev)
424	!c
425	INTEGER kmax, lmax, lmax3
426	PARAMETER(kmax=8, lmax=8, lmax3=3)
427	INTEGER kmaxm1, lmaxm1
428	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
429	INTEGER iimx7, jjmx7, jjmp1x7
430	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1, &
431	&jjmp1x7=jjmp1*lmaxm1)
432	!c
433	!c output from ISCCP simulator
434	REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
435	REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
436	REAL totalcldarea(klon,napisccp)
437	REAL meanptop(klon,napisccp)
438	REAL meantaucld(klon,napisccp)
439	REAL boxtau(klon,ncolmx,napisccp)
440	REAL boxptop(klon,ncolmx,napisccp)
441	REAL zx_tmp_fi3d_bx(klon,ncolmx)
442	REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
443	!c
444	REAL cld_fi3d(klon,lmax3)
445	REAL cld_3d(iim,jjmp1,lmax3)
446	!c
447	INTEGER iw, iwmax
448	REAL wmin, pas_w
449	!c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
450	!IM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
451	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
452	REAL o500(klon)
453	!c
454
455	!c sorties ISCCP
456
457	integer nid_isccp
458	save nid_isccp
459	!$OMP THREADPRIVATE(nid_isccp)
460
461	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
462	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
463	SAVE zx_tau
464	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
465	SAVE zx_pc
466	!$OMP THREADPRIVATE(zx_tau,zx_pc)
467	!c cldtopres pression au sommet des nuages
468	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
469	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
470	DATA cldtopres3/440., 680., 1000./
471	SAVE cldtopres,cldtopres3
472	!$OMP THREADPRIVATE(cldtopres,cldtopres3)
473	!IM 051005 BEG
474	INTEGER komega, nhoriRD
475
476	!c taulev: numero du niveau de tau dans les sorties ISCCP
477	CHARACTER *4 taulev(kmaxm1)
478	!c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
479	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
480	CHARACTER *3 pclev(lmaxm1)
481	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
482	SAVE taulev,pclev
483	!$OMP THREADPRIVATE(taulev,pclev)
484	!c
485	!c cnameisccp
486	CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
487	!IM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
488	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3', &
489	& 'pc= 180-310hPa, tau< 0.3', &
490	& 'pc= 310-440hPa, tau< 0.3', &
491	& 'pc= 440-560hPa, tau< 0.3', &
492	& 'pc= 560-680hPa, tau< 0.3', &
493	& 'pc= 680-800hPa, tau< 0.3', &
494	& 'pc= 800-1000hPa, tau< 0.3', &
495	& 'pc= 50-180hPa, tau= 0.3-1.3', &
496	& 'pc= 180-310hPa, tau= 0.3-1.3', &
497	& 'pc= 310-440hPa, tau= 0.3-1.3', &
498	& 'pc= 440-560hPa, tau= 0.3-1.3', &
499	& 'pc= 560-680hPa, tau= 0.3-1.3', &
500	& 'pc= 680-800hPa, tau= 0.3-1.3', &
501	& 'pc= 800-1000hPa, tau= 0.3-1.3', &
502	& 'pc= 50-180hPa, tau= 1.3-3.6', &
503	& 'pc= 180-310hPa, tau= 1.3-3.6', &
504	& 'pc= 310-440hPa, tau= 1.3-3.6', &
505	& 'pc= 440-560hPa, tau= 1.3-3.6', &
506	& 'pc= 560-680hPa, tau= 1.3-3.6', &
507	& 'pc= 680-800hPa, tau= 1.3-3.6', &
508	& 'pc= 800-1000hPa, tau= 1.3-3.6', &
509	& 'pc= 50-180hPa, tau= 3.6-9.4', &
510	& 'pc= 180-310hPa, tau= 3.6-9.4', &
511	& 'pc= 310-440hPa, tau= 3.6-9.4', &
512	& 'pc= 440-560hPa, tau= 3.6-9.4', &
513	& 'pc= 560-680hPa, tau= 3.6-9.4', &
514	& 'pc= 680-800hPa, tau= 3.6-9.4', &
515	& 'pc= 800-1000hPa, tau= 3.6-9.4', &
516	& 'pc= 50-180hPa, tau= 9.4-23', &
517	& 'pc= 180-310hPa, tau= 9.4-23', &
518	& 'pc= 310-440hPa, tau= 9.4-23', &
519	& 'pc= 440-560hPa, tau= 9.4-23', &
520	& 'pc= 560-680hPa, tau= 9.4-23', &
521	& 'pc= 680-800hPa, tau= 9.4-23', &
522	& 'pc= 800-1000hPa, tau= 9.4-23', &
523	& 'pc= 50-180hPa, tau= 23-60', &
524	& 'pc= 180-310hPa, tau= 23-60', &
525	& 'pc= 310-440hPa, tau= 23-60', &
526	& 'pc= 440-560hPa, tau= 23-60', &
527	& 'pc= 560-680hPa, tau= 23-60', &
528	& 'pc= 680-800hPa, tau= 23-60', &
529	& 'pc= 800-1000hPa, tau= 23-60', &
530	& 'pc= 50-180hPa, tau> 60.', &
531	& 'pc= 180-310hPa, tau> 60.', &
532	& 'pc= 310-440hPa, tau> 60.', &
533	& 'pc= 440-560hPa, tau> 60.', &
534	& 'pc= 560-680hPa, tau> 60.', &
535	& 'pc= 680-800hPa, tau> 60.', &
536	& 'pc= 800-1000hPa, tau> 60.'/
537	SAVE cnameisccp
538	!$OMP THREADPRIVATE(cnameisccp)
539	!c
540	!c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
541	!c INTEGER nhorix7
542	!IM: region='3d' <==> sorties en global
543	CHARACTER*3 region
544	PARAMETER(region='3d')
545	!c
546	!IM ISCCP simulator v3.4
547	!c
548	logical ok_hf
549	!c
550	integer nid_hf, nid_hf3d
551	save ok_hf, nid_hf, nid_hf3d
552	!$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
553	!c QUESTION : noms de variables ?
554
555	INTEGER longcles
556	PARAMETER ( longcles = 20 )
557	REAL clesphy0( longcles )
558	!c
559	!c Variables propres a la physique
560	INTEGER itap
561	SAVE itap ! compteur pour la physique
562	!$OMP THREADPRIVATE(itap)
563	!c
564	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
565	! real slp(klon) ! sea level pressure
566	!c
567	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
568	! REAL fevap(klon,nbsrf)
569	! REAL fluxlat(klon,nbsrf)
570	!c
571	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
572	! REAL qsol(klon)
573	REAL,save :: solarlong0
574	!$OMP THREADPRIVATE(solarlong0)
575
576	!c
577	!c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
578	!c
579	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
580	REAL zulow(klon),zvlow(klon)
581	!c
582	INTEGER igwd,idx(klon),itest(klon)
583	!c
584	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
585	! REAL agesno(klon,nbsrf)
586
587	!c
588	!c REAL,allocatable,save :: run_off_lic_0(:)
589	!cc$OMP THREADPRIVATE(run_off_lic_0)
590	!cym SAVE run_off_lic_0
591	!cKE43
592	!c Variables liees a la convection de K. Emanuel (sb):
593	!c
594	REAL bas, top ! cloud base and top levels
595	SAVE bas
596	SAVE top
597	!$OMP THREADPRIVATE(bas, top)
598
599	REAL wdn(klon), tdn(klon), qdn(klon)
600	!c
601	!c=================================================================================================
602	!cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
603	!c Variables li\'ees \`a la poche froide (jyg)
604
605	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
606	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
607	! REAL Vprecip(klon,klev+1) ! precipitation vertical profile
608	!c
609	REAL wape_prescr, fip_prescr
610	INTEGER it_wape_prescr
611	SAVE wape_prescr, fip_prescr, it_wape_prescr
612	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
613	!c
614	!c variables supplementaires de concvl
615	REAL Tconv(klon,klev)
616	REAL ment(klon,klev,klev),sij(klon,klev,klev)
617	REAL dd_t(klon,klev),dd_q(klon,klev)
618
619	real, save :: alp_bl_prescr=0.
620	real, save :: ale_bl_prescr=0.
621
622	real, save :: ale_max=1000.
623	real, save :: alp_max=2.
624
625	real, save :: wake_s_min_lsp=0.1
626
627	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
628	!$OMP THREADPRIVATE(ale_max,alp_max)
629	!$OMP THREADPRIVATE(wake_s_min_lsp)
630
631	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
632	! real ale_wake(klon)
633	! real alp_wake(klon)
634
635	real ok_wk_lsp(klon)
636
637	!cRC
638	!c Variables li\'ees \`a la poche froide (jyg et rr)
639	!c Version diagnostique pour l'instant : pas de r\'etroaction sur la convection
640
641	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
642
643	REAL wake_dth(klon,klev) ! wake : temp pot difference
644
645	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
646	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
647	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
648	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
649	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
650	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
651	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
652	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
653	! REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
654	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
655	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
656	REAL wake_spread(klon,klev) ! spreading term in wake_delt
657	!c
658	!cpourquoi y'a pas de save??
659	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
660	! REAL wake_h(klon) ! Wake : hauteur de la poche froide
661	!c
662	INTEGER wake_k(klon) ! Wake sommet
663	!c
664	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
665	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
666	!c
667	!cjyg
668	!ccc REAL wake_pe(klon) ! Wake potential energy - WAPE
669
670	REAL wake_gfl(klon) ! Gust Front Length
671	REAL wake_dens(klon)
672	!c
673	!c
674	REAL dt_dwn(klon,klev)
675	REAL dq_dwn(klon,klev)
676	REAL wdt_PBL(klon,klev)
677	REAL udt_PBL(klon,klev)
678	REAL wdq_PBL(klon,klev)
679	REAL udq_PBL(klon,klev)
680	REAL M_dwn(klon,klev)
681	REAL M_up(klon,klev)
682	REAL dt_a(klon,klev)
683	REAL dq_a(klon,klev)
684	REAL, SAVE :: alp_offset
685	!$OMP THREADPRIVATE(alp_offset)
686
687	!c
688	!cRR:fin declarations poches froides
689	!c=======================================================================================================
690
691	REAL zw2(klon,klev+1)
692
693	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
694	! REAL fraca(klon,klev+1)
695	REAL ztv(klon,klev),ztva(klon,klev)
696	REAL zpspsk(klon,klev)
697	REAL ztla(klon,klev),zqla(klon,klev)
698	REAL zthl(klon,klev)
699
700	!ccc nrlmd le 10/04/2012
701
702	!c--------Stochastic Boundary Layer Triggering: ALE_BL--------
703	!c---Propri\'et\'es du thermiques au LCL
704	real zlcl_th(klon) ! Altitude du LCL calcul\'e continument (pcon dans thermcell_main.F90)
705	real fraca0(klon) ! Fraction des thermiques au LCL
706	real w0(klon) ! Vitesse des thermiques au LCL
707	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
708	real tke0(klon,klev+1) ! TKE au dÃ©but du pas de temps
709	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
710	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
711
712	!c---Spectre de thermiques de type 2 au LCL
713	real n2(klon),s2(klon)
714	real ale_bl_stat(klon)
715
716	!c---D\'eclenchement stochastique
717	integer :: tau_trig(klon)
718	real proba_notrig(klon)
719	real random_notrig(klon)
720
721	!c--------Statistical Boundary Layer Closure: ALP_BL--------
722	!c---Profils de TKE dans et hors du thermique
723	real pbl_tke_input(klon,klev+1,nbsrf)
724	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
725	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
726
727	!c---Fermeture statistique
728	real alp_bl_det(klon) ! ALP d\'terministe du thermique unique
729	real alp_bl_fluct_m(klon) ! ALP li\'ee aux fluctuations de flux de masse sous-nuageux
730	real alp_bl_fluct_tke(klon) ! ALP li\'ee aux fluctuations d'\'energie cin\'etique sous-nuageuse
731	real alp_bl_conv(klon) ! ALP li\'ee \`a grande \'echelle
732	real alp_bl_stat(klon) ! ALP totale
733
734	!ccc fin nrlmd le 10/04/2012
735
736	!c Variables locales pour la couche limite (al1):
737	!c
738	!cAl1 REAL pblh(klon) ! Hauteur de couche limite
739	!cAl1 SAVE pblh
740	!c34EK
741	!c
742	!c Variables locales:
743	!c
744	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
745	! REAL cdragh(klon) ! drag coefficient pour T and Q
746	! REAL cdragm(klon) ! drag coefficient pour vent
747	!cAA
748	!cAA Pour phytrac
749	REAL u1(klon) ! vents dans la premiere couche U
750	REAL v1(klon) ! vents dans la premiere couche V
751
752	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule\n
753	! REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
754
755	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
756	!@$$ PARAMETER (offline=.false.)
757	!@$$ INTEGER physid
758	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
759	REAL frac_nucl(klon,klev) ! idem (nucleation)
760	! RomP >>>
761	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
762	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
763	! REAL beta_prec(klon,klev) ! taux de conv de l'eau cond (utilise)
764	! RomP <<<
765	INTEGER :: iii
766	REAL :: calday
767
768	!IM cf FH pour Tiedtke 080604
769	REAL rain_tiedtke(klon),snow_tiedtke(klon)
770	!c
771	!IM 050204 END
772	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
773	! REAL evap(klon), devap(klon) ! evaporation et sa derivee
774	REAL devap(klon) ! evaporation et sa derivee
775
776	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
777	! REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
778	REAL dsens(klon) ! chaleur sensible et sa derivee
779
780	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
781	! REAL bils(klon) ! bilan de chaleur au sol
782	! REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
783	!C ! type de sous-surface et pondere par la fraction
784	! REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
785	!C ! type de sous-surface et pondere par la fraction
786	REAL slab_wfbils(klon) ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
787
788	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
789	! REAL fder(klon)
790	! REAL ve(klon) ! integr. verticale du transport meri. de l'energie
791	! REAL vq(klon) ! integr. verticale du transport meri. de l'eau
792	! REAL ue(klon) ! integr. verticale du transport zonal de l'energie
793	! REAL uq(klon) ! integr. verticale du transport zonal de l'eau
794	!c
795
796	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
797	! REAL frugs(klon,nbsrf)
798	! REAL zxrugs(klon) ! longueur de rugosite
799	!c
800	!c Conditions aux limites
801	!c
802	!
803	REAL :: day_since_equinox
804	! Date de l'equinoxe de printemps
805	INTEGER, parameter :: mth_eq=3, day_eq=21
806	REAL :: jD_eq
807
808	LOGICAL, parameter :: new_orbit = .true.
809
810	!c
811	INTEGER lmt_pas
812	SAVE lmt_pas ! frequence de mise a jour
813	!$OMP THREADPRIVATE(lmt_pas)
814	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
815	! real zmasse(klon, llm),exner(klon, llm)
816	real exner(klon, llm)
817	!C (column-density of mass of air in a cell, in kg m-2)
818	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
819
820	!IM sorties
821	REAL un_jour
822	PARAMETER(un_jour=86400.)
823	!c======================================================================
824	!c
825	!c Declaration des procedures appelees
826	!c
827	EXTERNAL angle ! calculer angle zenithal du soleil
828	EXTERNAL alboc ! calculer l'albedo sur ocean
829	EXTERNAL ajsec ! ajustement sec
830	EXTERNAL conlmd ! convection (schema LMD)
831	!cKE43
832	EXTERNAL conema3 ! convect4.3
833	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
834	!cAA
835	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
836	!c ! stockage des coefficients necessaires au
837	!c ! lessivage OFF-LINE et ON-LINE
838	EXTERNAL hgardfou ! verifier les temperatures
839	EXTERNAL nuage ! calculer les proprietes radiatives
840	!CC EXTERNAL o3cm ! initialiser l'ozone
841	EXTERNAL orbite ! calculer l'orbite terrestre
842	EXTERNAL phyetat0 ! lire l'etat initial de la physique
843	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
844	EXTERNAL suphel ! initialiser certaines constantes
845	EXTERNAL transp ! transport total de l'eau et de l'energie
846	EXTERNAL ecribina ! ecrire le fichier binaire global
847	EXTERNAL ecribins ! ecrire le fichier binaire global
848	EXTERNAL ecrirega ! ecrire le fichier binaire regional
849	EXTERNAL ecriregs ! ecrire le fichier binaire regional
850	!IM
851	EXTERNAL haut2bas !variables de haut en bas
852	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
853	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
854	!c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
855	!c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
856	!c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
857	!c
858	!c Variables locales
859	!c
860	REAL rhcl(klon,klev) ! humiditi relative ciel clair
861	REAL dialiq(klon,klev) ! eau liquide nuageuse
862	REAL diafra(klon,klev) ! fraction nuageuse
863	REAL cldliq(klon,klev) ! eau liquide nuageuse
864	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
865	! REAL cldfra(klon,klev) ! fraction nuageuse
866	! REAL cldtau(klon,klev) ! epaisseur optique
867	! REAL cldemi(klon,klev) ! emissivite infrarouge
868	!c
869	!CXXX PB
870	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
871	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
872	! REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
873	! REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
874	! REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
875	!c
876	REAL zxfluxt(klon, klev)
877	REAL zxfluxq(klon, klev)
878	REAL zxfluxu(klon, klev)
879	REAL zxfluxv(klon, klev)
880	!CXXX
881	!c
882	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
883	! REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
884	! REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
885	!c Le rayonnement n'est pas calcule tous les pas, il faut donc
886	!c sauvegarder les sorties du rayonnement
887	!cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
888	!cym SAVE sollwdownclr, toplwdown, toplwdownclr
889	!cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
890	!c
891	INTEGER itaprad
892	SAVE itaprad
893	!$OMP THREADPRIVATE(itaprad)
894	!c
895	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
896	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
897	!c
898	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
899	! REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
900	! REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
901	!c
902	! REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
903	REAL zxsnow_dummy(klon)
904	!c
905	REAL dist, rmu0(klon), fract(klon)
906	REAL zdtime, zlongi
907	!c
908	CHARACTER*2 str2
909	CHARACTER*2 iqn
910	!c
911	REAL qcheck
912	REAL z_avant(klon), z_apres(klon), z_factor(klon)
913	LOGICAL zx_ajustq
914	!c
915	REAL za, zb
916	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
917	real zqsat(klon,klev)
918	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
919	REAL t_coup
920	PARAMETER (t_coup=234.0)
921	!c
922	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
923	! REAL zphi(klon,klev)
924	!cym A voir plus tard !!
925	!cym REAL zx_relief(iim,jjmp1)
926	!cym REAL zx_aire(iim,jjmp1)
927	!c
928	!c Grandeurs de sorties
929	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
930	! REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
931	! REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
932	! REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
933	REAL s_capCL(klon)
934	REAL s_oliqCL(klon), s_cteiCL(klon)
935	REAL s_trmb1(klon), s_trmb2(klon)
936	REAL s_trmb3(klon)
937	!cKE43
938	!c Variables locales pour la convection de K. Emanuel (sb):
939	!c
940	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
941	! REAL upwd(klon,klev) ! saturated updraft mass flux
942	! REAL dnwd(klon,klev) ! saturated downdraft mass flux
943	! REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
944	REAL tvp(klon,klev) ! virtual temp of lifted parcel
945	! REAL plcl(klon) ! Lifting Condensation Level
946	! REAL plfc(klon) ! Level of Free Convection
947	! REAL wbeff(klon) ! saturated updraft velocity at LFC
948	CHARACTER*40 capemaxcels !max(CAPE)
949
950	REAL rflag(klon) ! flag fonctionnement de convect
951	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
952
953	!c -- convect43:
954	INTEGER ntra ! nb traceurs pour convect4.3
955	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
956	REAL dtma_con(klon),dtlcl_con(klon)
957	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
958	REAL dplcldt(klon), dplcldr(klon)
959	!c? . condm_con(klon,klev),conda_con(klon,klev),
960	!c? . mr_con(klon,klev),ep_con(klon,klev)
961	!c? . ,sadiab(klon,klev),wadiab(klon,klev)
962	!c --
963	!c34EK
964	!c
965	!c Variables du changement
966	!c
967	!c con: convection
968	!c lsc: condensation a grande echelle (Large-Scale-Condensation)
969	!c ajs: ajustement sec
970	!c eva: evaporation de l'eau liquide nuageuse
971	!c vdf: couche limite (Vertical DiFfusion)
972	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
973	! REAL rneb(klon,klev)
974
975	! tendance nulles
976	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
977
978	!c
979	!*********************************************************
980	!* declarations
981
982	!*********************************************************
983	!IM 081204 END
984	!c
985	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
986	! REAL pmfu(klon,klev), pmfd(klon,klev)
987	REAL pen_u(klon,klev), pen_d(klon,klev)
988	REAL pde_u(klon,klev), pde_d(klon,klev)
989	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
990	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
991	! REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
992	! REAL prfl(klon,klev+1), psfl(klon,klev+1)
993	!c
994	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
995	! REAL rain_lsc(klon)
996	! REAL snow_lsc(klon)
997	!c
998	REAL ratqsc(klon,klev)
999	real ratqsbas,ratqshaut,tau_ratqs
1000	save ratqsbas,ratqshaut,tau_ratqs
1001	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
1002	real zpt_conv(klon,klev)
1003
1004	!c Parametres lies au nouveau schema de nuages (SB, PDF)
1005	real fact_cldcon
1006	real facttemps
1007	logical ok_newmicro
1008	save ok_newmicro
1009	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1010	! real ref_liq(klon,klev), ref_ice(klon,klev)
1011
1012	!$OMP THREADPRIVATE(ok_newmicro)
1013	save fact_cldcon,facttemps
1014	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
1015
1016	integer iflag_cldcon
1017	save iflag_cldcon
1018	!$OMP THREADPRIVATE(iflag_cldcon)
1019	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1020	! logical ptconv(klon,klev)
1021	!IM cf. AM 081204 BEG
1022	logical ptconvth(klon,klev)
1023	!IM cf. AM 081204 END
1024	!c
1025	!c Variables liees a l'ecriture de la bande histoire physique
1026	!c
1027	!c======================================================================
1028	!c
1029	!IM cf. AM 081204 BEG
1030	!c declarations pour sortir sur une sous-region
1031	integer imin_ins,imax_ins,jmin_ins,jmax_ins
1032	save imin_ins,imax_ins,jmin_ins,jmax_ins
1033	!$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
1034	!c real lonmin_ins,lonmax_ins,latmin_ins
1035	!c s ,latmax_ins
1036	!c data lonmin_ins,lonmax_ins,latmin_ins
1037	!c s ,latmax_ins/
1038	!c valeurs initiales s -5.,20.,41.,55./
1039	!c s 100.,130.,-20.,20./
1040	!c s -180.,180.,-90.,90./
1041	!c======================================================================
1042	!IM cf. AM 081204 END
1043
1044	!c
1045	integer itau_w ! pas de temps ecriture = itap + itau_phy
1046	!c
1047	!c
1048	!c Variables locales pour effectuer les appels en serie
1049	!c
1050	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule\n
1051	! REAL zx_rh(klon,klev)
1052	!IM RH a 2m (la surface)
1053	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1054	! REAL rh2m(klon), qsat2m(klon)
1055	! REAL tpot(klon), tpote(klon)
1056	REAL Lheat
1057
1058	INTEGER length
1059	PARAMETER ( length = 100 )
1060	REAL tabcntr0( length )
1061	!c
1062	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
1063	!IM
1064	INTEGER ndex2d1(iwmax)
1065	!c
1066	!IM AMIP2 BEG
1067	REAL moyglo, mountor
1068	!IM 141004 BEG
1069	REAL zustrdr(klon), zvstrdr(klon)
1070	REAL zustrli(klon), zvstrli(klon)
1071	REAL zustrph(klon), zvstrph(klon)
1072	REAL zustrhi(klon), zvstrhi(klon)
1073	REAL aam, torsfc
1074	!IM 141004 END
1075	!IM 190504 BEG
1076	INTEGER ij, imp1jmp1
1077	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
1078	!cym A voir plus tard
1079	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
1080	REAL padyn(iim+1,jjmp1,klev+1)
1081	REAL dudyn(iim+1,jjmp1,klev)
1082	REAL rlatdyn(iim+1,jjmp1)
1083	!IM 190504 END
1084	LOGICAL ok_msk
1085	REAL msk(klon)
1086	!IM
1087	REAL airetot, pi
1088	!cym A voir plus tard
1089	!cym REAL zm_wo(jjmp1, klev)
1090	!IM AMIP2 END
1091	!c
1092	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
1093	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
1094	REAL zx_tmp_fi3d1(klon,klev+1) !variable temporaire pour champs 3D (kelvp1)
1095	REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
1096	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
1097	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
1098	!c
1099	INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1100	INTEGER nid_hfnmc, nid_day_seri, nid_ctesGCM
1101	SAVE nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1102	SAVE nid_hfnmc, nid_day_seri, nid_ctesGCM
1103	!$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
1104	!$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
1105	!$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
1106	!c
1107	!IM 280405 BEG
1108	INTEGER nid_bilKPins, nid_bilKPave
1109	SAVE nid_bilKPins, nid_bilKPave
1110	!$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1111	!c
1112	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1113	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1114	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1115	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1116	!c
1117	INTEGER nhori, nvert, nvert1, nvert3
1118	REAL zsto, zsto1, zsto2
1119	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1120	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1121	REAL zout_isccp(napisccp)
1122	SAVE zcals, zcalh, zoutj, zout_isccp
1123	!$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1124
1125	real zjulian
1126	save zjulian
1127	!$OMP THREADPRIVATE(zjulian)
1128
1129	character*20 modname
1130	character*80 abort_message
1131	logical ok_sync
1132	real date0
1133	integer idayref
1134
1135	!C essai writephys
1136	integer fid_day, fid_mth, fid_ins
1137	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1138	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1139	parameter (prof2d_on = 1, prof3d_on = 2, &
1140	& prof2d_av = 3, prof3d_av = 4)
1141	character*30 nom_fichier
1142	character*40 varname
1143	character*40 vartitle
1144	character*20 varunits
1145	!C Variables liees au bilan d'energie et d'enthalpi
1146	REAL ztsol(klon)
1147	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
1148	& , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1149	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
1150	& , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1151	!$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1152	!$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1153	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1154	REAL d_h_vcol_phy
1155	REAL fs_bound, fq_bound
1156	SAVE d_h_vcol_phy
1157	!$OMP THREADPRIVATE(d_h_vcol_phy)
1158	REAL zero_v(klon)
1159	CHARACTER*15 ztit
1160	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1161	SAVE ip_ebil
1162	DATA ip_ebil/0/
1163	!$OMP THREADPRIVATE(ip_ebil)
1164	INTEGER if_ebil ! level for energy conserv. dignostics
1165	SAVE if_ebil
1166	!$OMP THREADPRIVATE(if_ebil)
1167	!c+jld ec_conser
1168	REAL ZRCPD
1169	!c-jld ec_conser
1170	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1171	! REAL t2m(klon,nbsrf) ! temperature a 2m
1172	REAL q2m(klon,nbsrf) ! humidite a 2m
1173
1174	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1175	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1176	! REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1177	! REAL zustar(klon),zu10m(klon), zv10m(klon) ! u* et vents a 10m moyennes s/1 maille
1178	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1179	CHARACTER*40 tinst, tave, typeval
1180	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1181
1182	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1183	! REAL re(klon, klev) ! Cloud droplet effective radius
1184	! REAL fl(klon, klev) ! denominator of re
1185
1186	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1187
1188	! Aerosol optical properties
1189	CHARACTER*4, DIMENSION(naero_grp) :: rfname
1190	REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index
1191	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass concentration for all soluble aerosols[ug/m3]
1192	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi ! - " - (pre-industrial value)
1193	INTEGER :: naero ! aerosol species
1194
1195	! Parameters
1196	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1197	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1198	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1199	SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
1200	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
1201	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1202	! false : lecture des aerosol dans un fichier
1203	!$OMP THREADPRIVATE(aerosol_couple)
1204	INTEGER, SAVE :: flag_aerosol
1205	!$OMP THREADPRIVATE(flag_aerosol)
1206	LOGICAL, SAVE :: new_aod
1207	!$OMP THREADPRIVATE(new_aod)
1208	!c
1209	!c--STRAT AEROSOL
1210	LOGICAL, SAVE :: flag_aerosol_strat
1211	!$OMP THREADPRIVATE(flag_aerosol_strat)
1212	!cc-fin STRAT AEROSOL
1213	!c
1214	!c Declaration des constantes et des fonctions thermodynamiques
1215	!c
1216	LOGICAL,SAVE :: first=.true.
1217	!$OMP THREADPRIVATE(first)
1218
1219	integer iunit
1220
1221	integer, save:: read_climoz ! read ozone climatology
1222	!C (let it keep the default OpenMP shared attribute)
1223	!C Allowed values are 0, 1 and 2
1224	!C 0: do not read an ozone climatology
1225	!C 1: read a single ozone climatology that will be used day and night
1226	!C 2: read two ozone climatologies, the average day and night
1227	!C climatology and the daylight climatology
1228
1229	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
1230	!C (let it keep the default OpenMP shared attribute)
1231
1232	real, pointer, save:: press_climoz(:)
1233	!C (let it keep the default OpenMP shared attribute)
1234	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
1235	! ascending order
1236
1237	integer, save:: co3i = 0
1238	! time index in NetCDF file of current ozone fields
1239	!$OMP THREADPRIVATE(co3i)
1240
1241	integer ro3i
1242	! required time index in NetCDF file for the ozone fields, between 1
1243	! and 360
1244
1245	INTEGER ierr
1246	#include "YOMCST.h"
1247	#include "YOETHF.h"
1248	#include "FCTTRE.h"
1249	!IM 100106 BEG : pouvoir sortir les ctes de la physique
1250	#include "conema3.h"
1251	#include "fisrtilp.h"
1252	#include "nuage.h"
1253	#include "compbl.h"
1254	!IM 100106 END : pouvoir sortir les ctes de la physique
1255	!c
1256	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1257	!c Declarations pour Simulateur COSP
1258	!c============================================================
1259	real :: mr_ozone(klon,klev)
1260
1261	!IM sorties fichier 1D paramLMDZ_phy.nc
1262	REAL :: zx_tmp_0d(1,1)
1263	INTEGER, PARAMETER :: np=1
1264	REAL,dimension(klon_glo) :: rlat_glo
1265	REAL,dimension(klon_glo) :: rlon_glo
1266	REAL gbils(1), gevap(1), gevapt(1), glat(1), gnet0(1), gnet(1)
1267	REAL grain(1), gtsol(1), gt2m(1), gprw(1)
1268
1269	!IM stations CFMIP
1270	INTEGER, SAVE :: nCFMIP
1271	!$OMP THREADPRIVATE(nCFMIP)
1272	INTEGER, PARAMETER :: npCFMIP=120
1273	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1274	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1275	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1276	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1277	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1278	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1279	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1280	!$OMP THREADPRIVATE(iGCM, jGCM)
1281	logical, dimension(nfiles) :: phys_out_filestations
1282	logical, parameter :: lNMC=.FALSE.
1283
1284	!IM betaCRF
1285	REAL, SAVE :: pfree, beta_pbl, beta_free
1286	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1287	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
1288	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
1289	LOGICAL, SAVE :: mskocean_beta
1290	!$OMP THREADPRIVATE(mskocean_beta)
1291	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
1292	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique pour radlwsw pour tester "CRF off"
1293	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique pour radlwsw pour tester "CRF off"
1294	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour radlwsw pour tester "CRF off"
1295	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
1296
1297	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1298	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1299	integer iostat
1300
1301	! L. Fita, LMD. Point for checkings...
1302	INTEGER :: llp
1303	llp = 644
1304
1305	! L. Fita, LMD August 2014
1306	CHARACTER(LEN=50) :: errmsg, fname, varname
1307	CHARACTER(LEN=50) :: jDS, jHS
1308	LOGICAL :: found
1309	REAL :: largest
1310
1311	fname = 'diagphy'
1312	errmsg = 'ERROR -- error -- ERROR -- error'
1313	largest = 10.e9
1314
1315	!c======================================================================
1316	! Gestion calendrier : mise a jour du module phys_cal_mod
1317	!
1318	CALL phys_cal_update(jD_cur,jH_cur)
1319	WRITE(jDS,'f10.5')jD_cur
1320	WRITE(jHS,'f10.5')jH_cur
1321
1322	fname = 'Entering in physic at day= '//TRIM(jDS)//' hour: '//TRIM(jHS)
1323	varname = 'paprs'
1324	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
1325	varname = 'pplay'
1326	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
1327	varname = 'pphi'
1328	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
1329	varname = 't_seri'
1330	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
1331	varname = 'u_seri'
1332	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
1333	varname = 'v_seri'
1334	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
1335	varname = 'q_seri'
1336	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
1337
1338
1339	!c======================================================================
1340	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1341	! Utilise notamment en 1D mais peut etre active egalement en 3D
1342	! en imposant la valeur de igout.
1343	!c======================================================================
1344
1345	if (prt_level.ge.1) then
1346	! Lluis
1347	! igout=klon/2+1/klon
1348	igout=llp
1349	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1350	write(lunout,*) &
1351	& 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1352	write(lunout,*) &
1353	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1354
1355	write(lunout,*) 'paprs, play, phi, u, v, t'
1356	do k=1,klev
1357	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1358	& u(igout,k),v(igout,k),t(igout,k)
1359	enddo
1360	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1361	do k=1,klev
1362	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1363	enddo
1364	endif
1365
1366	!c======================================================================
1367
1368	!cym => necessaire pour iflag_con != 2
1369	pmfd(:,:) = 0.
1370	pen_u(:,:) = 0.
1371	pen_d(:,:) = 0.
1372	pde_d(:,:) = 0.
1373	pde_u(:,:) = 0.
1374	aam=0.
1375
1376	torsfc=0.
1377
1378	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1379
1380	if (first) then
1381
1382	!cCR:nvelles variables convection/poches froides
1383
1384	print*, '================================================='
1385	print*, 'Allocation des variables locales et sauvegardees'
1386
1387	call phys_local_var_init
1388
1389	!c
1390	pasphys=pdtphys
1391	!c appel a la lecture du run.def physique
1392	call conf_phys(ok_journe, ok_mensuel, &
1393	& ok_instan, ok_hf, &
1394	& ok_LES, &
1395	& callstats, &
1396	& solarlong0,seuil_inversion, &
1397	& fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1398	& iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1399	& ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
1400	& flag_aerosol, flag_aerosol_strat, new_aod, &
1401	& bl95_b0, bl95_b1, &
1402	!c nv flags pour la convection et les poches froides &
1403	& read_climoz, &
1404	& alp_offset)
1405
1406
1407	!! call phys_state_var_init(read_climoz)
1408	call phys_output_var_init
1409
1410
1411	!L. Fita, LMD. November 2013
1412	! Initializing
1413	! Variables have been already initialilzed from the WRF variables.
1414	! Using on purpose subroutines from physiq_limit_variables_mod
1415	!! CALL neige_initialize()
1416	!! CALL limit_initialize()
1417	!! CALL vars_limit_init(klon)
1418
1419	! L. Fita, LMD January 2014. Initializing output variables
1420	CALL init_ovars_lmdz_NOmodule(klon,klev,nbsrf)
1421
1422	! Allocating restart variables
1423	! ALLOCATE(qsol(klon))
1424	! ALLOCATE(fder(klon))
1425	! ALLOCATE(snow(klon, nbsrf))
1426	! ALLOCATE(qsurf(klon, nbsrf))
1427	! ALLOCATE(evap(klon, nbsrf))
1428	! ALLOCATE(rugos(klon, nbsrf))
1429	! ALLOCATE(agesno(klon, nbsrf))
1430	! ALLOCATE(ftsoil(klon, nsoilmx, nbsrf))
1431	!! ALLOCATE(restart_runoff(klon))
1432
1433	! CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst, &
1434	! rugos_rst, agesno_rst, ftsoil_rst)
1435	! CALL fonte_neige_init(restart_runoff)
1436	! qsol=qsol_rst
1437
1438	print*, '================================================='
1439	!c
1440	dnwd0=0.0
1441	ftd=0.0
1442	fqd=0.0
1443	cin=0.
1444	!cym Attention pbase pas initialise dans concvl !!!!
1445	pbase=0
1446	!IM 180608
1447
1448	itau_con=0
1449	first=.false.
1450
1451	endif ! first
1452
1453	modname = 'physiq'
1454	!IM
1455	IF (ip_ebil_phy.ge.1) THEN
1456	DO i=1,klon
1457	zero_v(i)=0.
1458	END DO
1459	END IF
1460	ok_sync=.TRUE.
1461
1462	IF (debut) THEN
1463	CALL suphel ! initialiser constantes et parametres phys.
1464	ENDIF
1465
1466	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1467
1468
1469	!c======================================================================
1470	! Gestion calendrier : mise a jour du module phys_cal_mod
1471	!
1472	!c CALL phys_cal_update(jD_cur,jH_cur)
1473
1474	!c
1475	!c Si c'est le debut, il faut initialiser plusieurs choses
1476	!c ********
1477	!c
1478
1479	IF (debut) THEN
1480	!rv
1481	!cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1482	!cde la convection a partir des caracteristiques du thermique
1483	wght_th(:,:)=1.
1484	lalim_conv(:)=1
1485	!cRC
1486	ustar(:,:)=0.
1487	u10m(:,:)=0.
1488	v10m(:,:)=0.
1489	rain_con(:)=0.
1490	snow_con(:)=0.
1491	topswai(:)=0.
1492	topswad(:)=0.
1493	solswai(:)=0.
1494	solswad(:)=0.
1495
1496	wmax_th(:)=0.
1497	tau_overturning_th(:)=0.
1498
1499	IF (type_trac == 'inca') THEN
1500	! jg : initialisation jusqu'au ces variables sont dans restart
1501	ccm(:,:,:) = 0.
1502	tau_aero(:,:,:,:) = 0.
1503	piz_aero(:,:,:,:) = 0.
1504	cg_aero(:,:,:,:) = 0.
1505	END IF
1506
1507	rnebcon0(:,:) = 0.0
1508	clwcon0(:,:) = 0.0
1509	rnebcon(:,:) = 0.0
1510	clwcon(:,:) = 0.0
1511
1512	!IM
1513	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1514	!c
1515	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1516	& iflag_coupl,iflag_clos,iflag_wake
1517	print*,'CYCLE_DIURNE', cycle_diurne
1518	!c
1519	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1520	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1521	CALL abort_gcm (modname,abort_message,1)
1522	ENDIF
1523	!c
1524	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1525	abort_message = 'ISCCP-like outputs may be available for KE &
1526	&(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1527	CALL abort_gcm (modname,abort_message,1)
1528	ENDIF
1529	!c
1530	!c Initialiser les compteurs:
1531	!c
1532	itap = 0
1533	itaprad = 0
1534
1535	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1536	!! Un petit travail \`a faire ici.
1537	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1538
1539	if (iflag_pbl>1) then
1540	PRINT*, "Using method MELLOR&YAMADA"
1541	endif
1542
1543	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1544	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1545	! dyn3d
1546	! Attention : la version precedente n'etait pas tres propre.
1547	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1548	! pour obtenir le meme resultat.
1549	dtime=pdtphys
1550	radpas = NINT( 86400./dtime/nbapp_rad)
1551	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1552
1553	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1554
1555	IF (klon_glo==1) THEN
1556	coefh=0. ; coefm=0. ; pbl_tke=0.
1557	coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
1558	PRINT*,'FH WARNING : lignes a supprimer'
1559	ENDIF
1560	!IM begin
1561	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1562	&,ratqs(1,1)
1563	!IM end
1564
1565
1566
1567	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1568	!c
1569	!C on remet le calendrier a zero
1570	!c
1571	IF (raz_date .eq. 1) THEN
1572	itau_phy = 0
1573	ENDIF
1574
1575	!IM cf. AM 081204 BEG
1576	PRINT*,'cycle_diurne3 =',cycle_diurne
1577	!IM cf. AM 081204 END
1578	!c
1579	CALL printflag( tabcntr0,radpas,ok_journe, &
1580	& ok_instan, ok_region )
1581	!c
1582	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1583	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1584	& pdtphys
1585	abort_message='Pas physique n est pas correct '
1586	! call abort_gcm(modname,abort_message,1)
1587	dtime=pdtphys
1588	ENDIF
1589	IF (nlon .NE. klon) THEN
1590	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1591	& klon
1592	abort_message='nlon et klon ne sont pas coherents'
1593	call abort_gcm(modname,abort_message,1)
1594	ENDIF
1595	IF (nlev .NE. klev) THEN
1596	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1597	& klev
1598	abort_message='nlev et klev ne sont pas coherents'
1599	call abort_gcm(modname,abort_message,1)
1600	ENDIF
1601	!c
1602	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1603	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1604	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1605	abort_message='Nbre d appels au rayonnement insuffisant'
1606	call abort_gcm(modname,abort_message,1)
1607	ENDIF
1608	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1609	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1610	& ok_cvl
1611	!c
1612	!cKE43
1613	!c Initialisation pour la convection de K.E. (sb):
1614	IF (iflag_con.GE.3) THEN
1615
1616	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1617	WRITE(lunout,*) &
1618	& "On va utiliser le melange convectif des traceurs qui"
1619	WRITE(lunout,*)"est calcule dans convect4.3"
1620	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1621
1622	DO i = 1, klon
1623	ema_cbmf(i) = 0.
1624	ema_pcb(i) = 0.
1625	ema_pct(i) = 0.
1626	!c ema_workcbmf(i) = 0.
1627	ENDDO
1628	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1629	DO i = 1, klon
1630	ibas_con(i) = 1
1631	itop_con(i) = 1
1632	ENDDO
1633	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1634	!c===============================================================================
1635	!cCR:04.12.07: initialisations poches froides
1636	!c Controle de ALE et ALP pour la fermeture convective (jyg)
1637	if (iflag_wake>=1) then
1638	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1639	& ,alp_bl_prescr, ale_bl_prescr)
1640	!c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1641	!c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1642	endif
1643
1644	do i = 1,klon
1645	Ale_bl(i)=0.
1646	Alp_bl(i)=0.
1647	enddo
1648
1649	!c================================================================================
1650	!IM stations CFMIP
1651	nCFMIP=npCFMIP
1652	OPEN(98,file='npCFMIP_param.data',status='old', &
1653	& form='formatted',iostat=iostat)
1654	if (iostat == 0) then
1655	READ(98,*,end=998) nCFMIP
1656	998 CONTINUE
1657	CLOSE(98)
1658	CONTINUE
1659	IF(nCFMIP.GT.npCFMIP) THEN
1660	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1661	CALL abort
1662	else
1663	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1664	ENDIF
1665
1666	!c
1667	ALLOCATE(tabCFMIP(nCFMIP))
1668	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1669	ALLOCATE(tabijGCM(nCFMIP))
1670	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1671	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1672	!c
1673	!c lecture des nCFMIP stations CFMIP, de leur numero
1674	!c et des coordonnees geographiques lonCFMIP, latCFMIP
1675	!c
1676	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1677	&lonCFMIP, latCFMIP)
1678	!c
1679	!c identification des
1680	!c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
1681	!c 2) indices points tabijGCM de la grille physique 1d sur klon points
1682	!c 3) indices iGCM, jGCM de la grille physique 2d
1683	!c
1684	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1685	&tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1686	!c
1687	else
1688	ALLOCATE(tabijGCM(0))
1689	ALLOCATE(lonGCM(0), latGCM(0))
1690	ALLOCATE(iGCM(0), jGCM(0))
1691	end if
1692	else
1693	ALLOCATE(tabijGCM(0))
1694	ALLOCATE(lonGCM(0), latGCM(0))
1695	ALLOCATE(iGCM(0), jGCM(0))
1696	ENDIF
1697
1698	DO i=1,klon
1699	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1700	ENDDO
1701
1702	!c34EK
1703	IF (ok_orodr) THEN
1704
1705	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1706	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1707	! justement quand ok_orodr = false.
1708	! ce rugoro est utilise par la couche limite et fait double emploi
1709	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1710	! DO i=1,klon
1711	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1712	! ENDDO
1713	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1714	IF (ok_strato) THEN
1715	CALL SUGWD_strato(klon,klev,paprs,pplay)
1716	ELSE
1717	CALL SUGWD(klon,klev,paprs,pplay)
1718	ENDIF
1719
1720	DO i=1,klon
1721	zuthe(i)=0.
1722	zvthe(i)=0.
1723	if(zstd(i).gt.10.)then
1724	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1725	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1726	endif
1727	ENDDO
1728	ENDIF
1729	!c
1730	!c
1731	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1732	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1733	& lmt_pas
1734	!c
1735	capemaxcels = 't_max(X)'
1736	t2mincels = 't_min(X)'
1737	t2maxcels = 't_max(X)'
1738	tinst = 'inst(X)'
1739	tave = 'ave(X)'
1740	!IM cf. AM 081204 BEG
1741	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1742	!IM cf. AM 081204 END
1743	!c
1744	!c=============================================================
1745	!c Initialisation des sorties
1746	!c=============================================================
1747
1748	#ifdef CPP_IOIPSL
1749
1750	!$OMP MASTER
1751	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
1752	& iGCM,jGCM,lonGCM,latGCM, &
1753	& jjmp1,nlevSTD,clevSTD, &
1754	& nbteta, ctetaSTD, dtime,ok_veget, &
1755	& type_ocean,iflag_pbl,ok_mensuel,ok_journe, &
1756	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1757	& read_climoz, phys_out_filestations, &
1758	& new_aod, aerosol_couple, &
1759	& flag_aerosol_strat )
1760	!$OMP END MASTER
1761	!$OMP BARRIER
1762
1763	#ifdef histISCCP
1764	#include "ini_histISCCP.h"
1765	#endif
1766
1767	#ifdef histNMC
1768	#include "ini_histhfNMC.h"
1769	#include "ini_histdayNMC.h"
1770	#include "ini_histmthNMC.h"
1771	#endif
1772
1773	#include "ini_histday_seri.h"
1774
1775	#include "ini_paramLMDZ_phy.h"
1776
1777	#endif
1778	ecrit_reg = ecrit_reg * un_jour
1779	ecrit_tra = ecrit_tra * un_jour
1780
1781	!cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1782	date0 = jD_ref
1783	WRITE(,) 'physiq date0 : ',date0
1784	!c
1785	!c
1786	!c
1787	!c Prescrire l'ozone dans l'atmosphere
1788	!c
1789	!c
1790	!cc DO i = 1, klon
1791	!cc DO k = 1, klev
1792	!cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1793	!cc ENDDO
1794	!cc ENDDO
1795	!c
1796	IF (type_trac == 'inca') THEN
1797	#ifdef INCA
1798	CALL VTe(VTphysiq)
1799	CALL VTb(VTinca)
1800	! iii = MOD(NINT(xjour),360)
1801	! calday = REAL(iii) + jH_cur
1802	calday = REAL(days_elapsed) + jH_cur
1803	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1804
1805	CALL chemini( &
1806	& rg, &
1807	& ra, &
1808	& airephy, &
1809	& rlat, &
1810	& rlon, &
1811	& presnivs, &
1812	& calday, &
1813	& klon, &
1814	& nqtot, &
1815	& pdtphys, &
1816	& annee_ref, &
1817	& day_ref, &
1818	& itau_phy)
1819
1820	CALL VTe(VTinca)
1821	CALL VTb(VTphysiq)
1822	#endif
1823	END IF
1824	!c
1825	!c
1826	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1827	! Nouvelle initialisation pour le rayonnement RRTM
1828	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1829
1830	call iniradia(klon,klev,paprs(1,1:klev+1))
1831
1832	!$OMP single
1833	if (read_climoz >= 1) then
1834	call open_climoz(ncid_climoz, press_climoz)
1835	END IF
1836	!$OMP end single
1837	!c
1838	!IM betaCRF
1839	pfree=70000. !Pa
1840	beta_pbl=1.
1841	beta_free=1.
1842	lon1_beta=-180.
1843	lon2_beta=+180.
1844	lat1_beta=90.
1845	lat2_beta=-90.
1846	mskocean_beta=.FALSE.
1847
1848	OPEN(99,file='beta_crf.data',status='old', &
1849	& form='formatted',err=9999)
1850	READ(99,*,end=9998) pfree
1851	READ(99,*,end=9998) beta_pbl
1852	READ(99,*,end=9998) beta_free
1853	READ(99,*,end=9998) lon1_beta
1854	READ(99,*,end=9998) lon2_beta
1855	READ(99,*,end=9998) lat1_beta
1856	READ(99,*,end=9998) lat2_beta
1857	READ(99,*,end=9998) mskocean_beta
1858	9998 Continue
1859	CLOSE(99)
1860	9999 Continue
1861	WRITE(,)'pfree=',pfree
1862	WRITE(,)'beta_pbl=',beta_pbl
1863	WRITE(,)'beta_free=',beta_free
1864	WRITE(,)'lon1_beta=',lon1_beta
1865	WRITE(,)'lon2_beta=',lon2_beta
1866	WRITE(,)'lat1_beta=',lat1_beta
1867	WRITE(,)'lat2_beta=',lat2_beta
1868	WRITE(,)'mskocean_beta=',mskocean_beta
1869	ENDIF
1870	!
1871	! ************** Fin de IF ( debut ) *************
1872	!
1873	!
1874	! Incrementer le compteur de la physique
1875	!
1876	itap = itap + 1
1877	WRITE(jDS,'I5')itap
1878	fname = 'After itap= '//TRIM(jDS)
1879	varname = 'paprs'
1880	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
1881	varname = 'pplay'
1882	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
1883	varname = 'pphi'
1884	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
1885	varname = 't_seri'
1886	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
1887	varname = 'u_seri'
1888	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
1889	varname = 'v_seri'
1890	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
1891	varname = 'q_seri'
1892	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
1893
1894	!c
1895	!
1896	! Update fraction of the sub-surfaces (pctsrf) and
1897	! initialize, where a new fraction has appeared, all variables depending
1898	! on the surface fraction.
1899	!
1900	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1901	& pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1902
1903	! Update time and other variables in Reprobus
1904	IF (type_trac == 'repr') THEN
1905	#ifdef REPROBUS
1906	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1907	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1908	CALL Rtime(debut)
1909	#endif
1910	END IF
1911
1912
1913	! Tendances bidons pour les processus qui n'affectent pas certaines
1914	! variables.
1915	du0(:,:)=0.
1916	dv0(:,:)=0.
1917	dq0(:,:)=0.
1918	dql0(:,:)=0.
1919	!c
1920	!c Mettre a zero des variables de sortie (pour securite)
1921	!c
1922	DO i = 1, klon
1923	d_ps(i) = 0.0
1924	ENDDO
1925	DO k = 1, klev
1926	DO i = 1, klon
1927	d_t(i,k) = 0.0
1928	d_u(i,k) = 0.0
1929	d_v(i,k) = 0.0
1930	ENDDO
1931	ENDDO
1932	DO iq = 1, nqtot
1933	DO k = 1, klev
1934	DO i = 1, klon
1935	d_qx(i,k,iq) = 0.0
1936	ENDDO
1937	ENDDO
1938	ENDDO
1939	da(:,:)=0.
1940	mp(:,:)=0.
1941	phi(:,:,:)=0.
1942	! RomP >>>
1943	phi2(:,:,:)=0.
1944	beta_prec_fisrt(:,:)=0.
1945	beta_prec(:,:)=0.
1946	epmlmMm(:,:,:)=0.
1947	eplaMm(:,:)=0.
1948	d1a(:,:)=0.
1949	dam(:,:)=0.
1950	pmflxr=0.
1951	pmflxs=0.
1952	! RomP <<<
1953
1954	!c
1955	!c Ne pas affecter les valeurs entrees de u, v, h, et q
1956	!c
1957	DO k = 1, klev
1958	DO i = 1, klon
1959	t_seri(i,k) = t(i,k)
1960	u_seri(i,k) = u(i,k)
1961	v_seri(i,k) = v(i,k)
1962	q_seri(i,k) = qx(i,k,ivap)
1963	ql_seri(i,k) = qx(i,k,iliq)
1964	qs_seri(i,k) = 0.
1965	ENDDO
1966	ENDDO
1967	tke0(:,:)=pbl_tke(:,:,is_ave)
1968	IF (nqtot.GE.3) THEN
1969	DO iq = 3, nqtot
1970	DO k = 1, klev
1971	DO i = 1, klon
1972	tr_seri(i,k,iq-2) = qx(i,k,iq)
1973	ENDDO
1974	ENDDO
1975	ENDDO
1976	ELSE
1977	DO k = 1, klev
1978	DO i = 1, klon
1979	tr_seri(i,k,1) = 0.0
1980	ENDDO
1981	ENDDO
1982	ENDIF
1983	!C
1984	DO i = 1, klon
1985	ztsol(i) = 0.
1986	ENDDO
1987	DO nsrf = 1, nbsrf
1988	DO i = 1, klon
1989	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1990	ENDDO
1991	ENDDO
1992	!IM
1993	IF (ip_ebil_phy.ge.1) THEN
1994	ztit='after dynamic'
1995	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
1996	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1997	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1998	!C Comme les tendances de la physique sont ajoute dans la dynamique,
1999	!C on devrait avoir que la variation d'entalpie par la dynamique
2000	!C est egale a la variation de la physique au pas de temps precedent.
2001	!C Donc la somme de ces 2 variations devrait etre nulle.
2002	call diagphy(airephy,ztit,ip_ebil_phy &
2003	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2004	& , zero_v, zero_v, zero_v, ztsol &
2005	& , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
2006	& , fs_bound, fq_bound )
2007	END IF
2008	PRINT *,' Lluis 1 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2009	fname = 'After dynamic'
2010	varname = 'paprs'
2011	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
2012	varname = 'pplay'
2013	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
2014	varname = 'pphi'
2015	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
2016	varname = 't_seri'
2017	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
2018	varname = 'u_seri'
2019	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
2020	varname = 'v_seri'
2021	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
2022	varname = 'q_seri'
2023	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
2024
2025	!c Diagnostiquer la tendance dynamique
2026	!c
2027	IF (ancien_ok) THEN
2028	DO k = 1, klev
2029	DO i = 1, klon
2030	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
2031	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
2032	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
2033	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
2034	ENDDO
2035	ENDDO
2036	!!! RomP >>> td dyn traceur
2037	IF (nqtot.GE.3) THEN
2038	DO iq = 3, nqtot
2039	DO k = 1, klev
2040	DO i = 1, klon
2041	d_tr_dyn(i,k,iq-2)= &
2042	& (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
2043	! iiq=niadv(iq)
2044	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
2045	ENDDO
2046	ENDDO
2047	ENDDO
2048	ENDIF
2049	!!! RomP <<<
2050	ELSE
2051	DO k = 1, klev
2052	DO i = 1, klon
2053	d_u_dyn(i,k) = 0.0
2054	d_v_dyn(i,k) = 0.0
2055	d_t_dyn(i,k) = 0.0
2056	d_q_dyn(i,k) = 0.0
2057	ENDDO
2058	ENDDO
2059	!!! RomP >>> td dyn traceur
2060	IF (nqtot.GE.3) THEN
2061	DO iq = 3, nqtot
2062	DO k = 1, klev
2063	DO i = 1, klon
2064	d_tr_dyn(i,k,iq-2)= 0.0
2065	ENDDO
2066	ENDDO
2067	ENDDO
2068	ENDIF
2069	!!! RomP <<<
2070	ancien_ok = .TRUE.
2071	ENDIF
2072	!c
2073	!c Ajouter le geopotentiel du sol:
2074	!c
2075	DO k = 1, klev
2076	DO i = 1, klon
2077	zphi(i,k) = pphi(i,k) + pphis(i)
2078	ENDDO
2079	ENDDO
2080	!c
2081	!c Verifier les temperatures
2082	!c
2083	!IM BEG
2084	IF (check) THEN
2085	amn=MIN(ftsol(1,is_ter),1000.)
2086	amx=MAX(ftsol(1,is_ter),-1000.)
2087	DO i=2, klon
2088	amn=MIN(ftsol(i,is_ter),amn)
2089	amx=MAX(ftsol(i,is_ter),amx)
2090	ENDDO
2091	!c
2092	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2093	ENDIF !(check) THEN
2094
2095	!IM END
2096	!c
2097
2098	CALL hgardfou(t_seri,ftsol,'debutphy')
2099	!c
2100	!IM BEG
2101	IF (check) THEN
2102	amn=MIN(ftsol(1,is_ter),1000.)
2103	amx=MAX(ftsol(1,is_ter),-1000.)
2104	DO i=2, klon
2105	amn=MIN(ftsol(i,is_ter),amn)
2106	amx=MAX(ftsol(i,is_ter),amx)
2107	ENDDO
2108	!c
2109	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2110	ENDIF !(check) THEN
2111
2112	!IM END
2113	!c
2114	!c Mettre en action les conditions aux limites (albedo, sst, etc.).
2115	!c Prescrire l'ozone et calculer l'albedo sur l'ocean.
2116	!c
2117	if (read_climoz >= 1) then
2118	!C Ozone from a file
2119	! Update required ozone index:
2120	ro3i = int((days_elapsed + jh_cur - jh_1jan) &
2121	& / ioget_year_len(year_cur) * 360.) + 1
2122	if (ro3i == 361) ro3i = 360
2123	!C (This should never occur, except perhaps because of roundup
2124	!C error. See documentation.)
2125	if (ro3i /= co3i) then
2126	!C Update ozone field:
2127	if (read_climoz == 1) then
2128	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
2129	& press_in_edg=press_climoz, paprs=paprs, v3=wo)
2130	else
2131	!C read_climoz == 2
2132	call regr_pr_av(ncid_climoz, &
2133	& (/"tro3 ", "tro3_daylight"/), &
2134	& julien=ro3i, press_in_edg=press_climoz, paprs=paprs, &
2135	& v3=wo)
2136	end if
2137	! Convert from mole fraction of ozone to column density of ozone in a
2138	! cell, in kDU:
2139	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) &
2140	& * rmo3 / rmd * zmasse / dobson_u / 1e3
2141	!C (By regridding ozone values for LMDZ only once every 360th of
2142	!C year, we have already neglected the variation of pressure in one
2143	!C 360th of year. So do not recompute "wo" at each time step even if
2144	!C "zmasse" changes a little.)
2145	co3i = ro3i
2146	end if
2147	elseif (MOD(itap-1,lmt_pas) == 0) THEN
2148	!C Once per day, update ozone from Royer:
2149	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
2150	ENDIF
2151
2152	!c
2153	!c Re-evaporer l'eau liquide nuageuse
2154	!c
2155	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
2156	DO i = 1, klon
2157	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2158	!c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2159	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2160	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
2161	zb = MAX(0.0,ql_seri(i,k))
2162	za = - MAX(0.0,ql_seri(i,k)) &
2163	& * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
2164	t_seri(i,k) = t_seri(i,k) + za
2165	q_seri(i,k) = q_seri(i,k) + zb
2166	ql_seri(i,k) = 0.0
2167	d_t_eva(i,k) = za
2168	d_q_eva(i,k) = zb
2169	ENDDO
2170	ENDDO
2171	!IM
2172	IF (ip_ebil_phy.ge.2) THEN
2173	ztit='after reevap'
2174	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
2175	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2176	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2177	call diagphy(airephy,ztit,ip_ebil_phy &
2178	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2179	& , zero_v, zero_v, zero_v, ztsol &
2180	& , d_h_vcol, d_qt, d_ec &
2181	& , fs_bound, fq_bound )
2182	!C
2183	END IF
2184	PRINT *,' Lluis 2 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2185	fname = 'After reevap'
2186	varname = 'paprs'
2187	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
2188	varname = 'pplay'
2189	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
2190	varname = 'pphi'
2191	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
2192	varname = 't_seri'
2193	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
2194	varname = 'u_seri'
2195	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
2196	varname = 'v_seri'
2197	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
2198	varname = 'q_seri'
2199	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
2200	!c
2201	!c=========================================================================
2202	! Calculs de l'orbite.
2203	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2204	! doit donc etre plac\'e avant radlwsw et pbl_surface
2205
2206	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2207	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2208	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2209	!
2210	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2211	! solarlong0
2212	if (solarlong0<-999.) then
2213	if (new_orbit) then
2214	! calcul selon la routine utilisee pour les planetes
2215	call solarlong(day_since_equinox, zlongi, dist)
2216	else
2217	! calcul selon la routine utilisee pour l'AR4
2218	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2219	endif
2220	else
2221	zlongi=solarlong0 ! longitude solaire vraie
2222	dist=1. ! distance au soleil / moyenne
2223	endif
2224	if(prt_level.ge.1) &
2225	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist,&
2226	' jD_cur: ',jD_cur,' jH_cur: ',jH_cur,' days_elapsed: ', &
2227	REAL(days_elapsed+1)
2228
2229
2230	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2231	! Calcul de l'ensoleillement :
2232	! ============================
2233	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2234	! l'annee a partir d'une formule analytique.
2235	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2236	! non nul aux poles.
2237	IF (abs(solarlong0-1000.)<1.e-4) then
2238	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
2239	ELSE
2240	! Avec ou sans cycle diurne
2241	IF (cycle_diurne) THEN
2242	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
2243	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
2244	ELSE
2245	CALL angle(zlongi, rlat, fract, rmu0)
2246	ENDIF
2247	ENDIF
2248
2249	if (mydebug) then
2250	call writefield_phy('u_seri',u_seri,llm)
2251	call writefield_phy('v_seri',v_seri,llm)
2252	call writefield_phy('t_seri',t_seri,llm)
2253	call writefield_phy('q_seri',q_seri,llm)
2254	endif
2255
2256	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2257	!c Appel au pbl_surface : Planetary Boudary Layer et Surface
2258	!c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2259	!c turbulent du couche limit.
2260	!c
2261	!c Certains varibales de sorties de pbl_surface sont utiliser que pour
2262	!c ecriture des fihiers hist_XXXX.nc, ces sont :
2263	!c qsol, zq2m, s_pblh, s_lcl,
2264	!c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2265	!c s_therm, s_trmb1, s_trmb2, s_trmb3,
2266	!c zxrugs, zu10m, zv10m, fder,
2267	!c zxqsurf, rh2m, zxfluxu, zxfluxv,
2268	!c frugs, agesno, fsollw, fsolsw,
2269	!c d_ts, fevap, fluxlat, t2m,
2270	!c wfbils, wfbilo, fluxt, fluxu, fluxv,
2271	!c
2272	!c Certains ne sont pas utiliser du tout :
2273	!c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2274	!c
2275
2276	!c Calcul de l'humidite de saturation au niveau du sol
2277
2278
2279
2280	if (iflag_pbl/=0) then
2281
2282	CALL pbl_surface( &
2283	& dtime, date0, itap, days_elapsed+1, &
2284	& debut, lafin, &
2285	& rlon, rlat, rugoro, rmu0, &
2286	& rain_fall, snow_fall, solsw, sollw, &
2287	& t_seri, q_seri, u_seri, v_seri, &
2288	& pplay, paprs, pctsrf, &
2289	& ftsol, falb1, falb2, ustar, u10m, v10m, &
2290	& sollwdown, cdragh, cdragm, u1, v1, &
2291	& albsol1, albsol2, sens, evap, &
2292	& zxtsol, zxfluxlat, zt2m, qsat2m, &
2293	& d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2294	& coefh, coefm, slab_wfbils, &
2295	& qsol, zq2m, s_pblh, s_lcl, &
2296	& s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2297	& s_therm, s_trmb1, s_trmb2, s_trmb3, &
2298	& zxrugs, zustar, zu10m, zv10m, fder, &
2299	& zxqsurf, rh2m, zxfluxu, zxfluxv, &
2300	& frugs, agesno, fsollw, fsolsw, &
2301	& d_ts, fevap, fluxlat, t2m, &
2302	& wfbils, wfbilo, fluxt, fluxu, fluxv, &
2303	& dsens, devap, zxsnow, &
2304	& zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
2305
2306	fname = 'After pbl_surface'
2307	varname = 'paprs'
2308	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
2309	varname = 'pplay'
2310	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
2311	varname = 'pphi'
2312	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
2313	varname = 't_seri'
2314	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
2315	varname = 'u_seri'
2316	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
2317	varname = 'v_seri'
2318	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
2319	varname = 'q_seri'
2320	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
2321
2322	!-----------------------------------------------------------------------------------------
2323	! ajout des tendances de la diffusion turbulente
2324	CALL add_phys_tend &
2325	& (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
2326	!-----------------------------------------------------------------------------------------
2327	if (mydebug) then
2328	call writefield_phy('u_seri',u_seri,llm)
2329	call writefield_phy('v_seri',v_seri,llm)
2330	call writefield_phy('t_seri',t_seri,llm)
2331	call writefield_phy('q_seri',q_seri,llm)
2332	endif
2333
2334	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2335	& t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2336
2337
2338	IF (ip_ebil_phy.ge.2) THEN
2339	ztit='after surface_main'
2340	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2341	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2342	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2343	call diagphy(airephy,ztit,ip_ebil_phy &
2344	& , zero_v, zero_v, zero_v, zero_v, sens &
2345	& , evap , zero_v, zero_v, ztsol &
2346	& , d_h_vcol, d_qt, d_ec &
2347	& , fs_bound, fq_bound )
2348	END IF
2349
2350	ENDIF
2351	PRINT *,' Lluis 3 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2352	fname = 'After surface_main'
2353	varname = 'paprs'
2354	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
2355	varname = 'pplay'
2356	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
2357	varname = 'pphi'
2358	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
2359	varname = 't_seri'
2360	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
2361	varname = 'u_seri'
2362	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
2363	varname = 'v_seri'
2364	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
2365	varname = 'q_seri'
2366	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
2367
2368	!c =================================================================== c
2369	!c Calcul de Qsat
2370
2371	DO k = 1, klev
2372	DO i = 1, klon
2373	zx_t = t_seri(i,k)
2374	IF (thermcep) THEN
2375	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2376	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2377	zx_qs = MIN(0.5,zx_qs)
2378	zcor = 1./(1.-retv*zx_qs)
2379	zx_qs = zx_qs*zcor
2380	ELSE
2381	IF (zx_t.LT.t_coup) THEN
2382	zx_qs = qsats(zx_t)/pplay(i,k)
2383	ELSE
2384	zx_qs = qsatl(zx_t)/pplay(i,k)
2385	ENDIF
2386	ENDIF
2387	zqsat(i,k)=zx_qs
2388	ENDDO
2389	ENDDO
2390
2391	if (prt_level.ge.1) then
2392	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2393	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2394	endif
2395	!c
2396	!c Appeler la convection (au choix)
2397	!c
2398	DO k = 1, klev
2399	DO i = 1, klon
2400	conv_q(i,k) = d_q_dyn(i,k) &
2401	& + d_q_vdf(i,k)/dtime
2402	conv_t(i,k) = d_t_dyn(i,k) &
2403	& + d_t_vdf(i,k)/dtime
2404	ENDDO
2405	ENDDO
2406	IF (check) THEN
2407	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2408	WRITE(lunout,*) "avantcon=", za
2409	ENDIF
2410	zx_ajustq = .FALSE.
2411	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2412	IF (zx_ajustq) THEN
2413	DO i = 1, klon
2414	z_avant(i) = 0.0
2415	ENDDO
2416	DO k = 1, klev
2417	DO i = 1, klon
2418	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2419	& *(paprs(i,k)-paprs(i,k+1))/RG
2420	ENDDO
2421	ENDDO
2422	ENDIF
2423
2424	!c Calcule de vitesse verticale a partir de flux de masse verticale
2425	DO k = 1, klev
2426	DO i = 1, klon
2427	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2428	END DO
2429	END DO
2430	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2431	& omega(igout, :)
2432
2433	IF (iflag_con.EQ.1) THEN
2434	abort_message ='reactiver le call conlmd dans physiq.F'
2435	CALL abort_gcm (modname,abort_message,1)
2436	!c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2437	!c . d_t_con, d_q_con,
2438	!c . rain_con, snow_con, ibas_con, itop_con)
2439	ELSE IF (iflag_con.EQ.2) THEN
2440	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
2441	& conv_t, conv_q, -evap, omega, &
2442	& d_t_con, d_q_con, rain_con, snow_con, &
2443	& pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2444	& kcbot, kctop, kdtop, pmflxr, pmflxs)
2445	d_u_con = 0.
2446	d_v_con = 0.
2447
2448	WHERE (rain_con < 0.) rain_con = 0.
2449	WHERE (snow_con < 0.) snow_con = 0.
2450	DO i = 1, klon
2451	ibas_con(i) = klev+1 - kcbot(i)
2452	itop_con(i) = klev+1 - kctop(i)
2453	ENDDO
2454	ELSE IF (iflag_con.GE.3) THEN
2455	!c nb of tracers for the KE convection:
2456	!c MAF la partie traceurs est faite dans phytrac
2457	!c on met ntra=1 pour limiter les appels mais on peut
2458	!c supprimer les calculs / ftra.
2459	ntra = 1
2460
2461	!c=====================================================================================
2462	!cajout pour la parametrisation des poches froides:
2463	!ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2464	do k=1,klev
2465	do i=1,klon
2466	if (iflag_wake>=1) then
2467	t_wake(i,k) = t_seri(i,k) &
2468	& +(1-wake_s(i))*wake_deltat(i,k)
2469	q_wake(i,k) = q_seri(i,k) &
2470	& +(1-wake_s(i))*wake_deltaq(i,k)
2471	t_undi(i,k) = t_seri(i,k) &
2472	& -wake_s(i)*wake_deltat(i,k)
2473	q_undi(i,k) = q_seri(i,k) &
2474	& -wake_s(i)*wake_deltaq(i,k)
2475	else
2476	t_wake(i,k) = t_seri(i,k)
2477	q_wake(i,k) = q_seri(i,k)
2478	t_undi(i,k) = t_seri(i,k)
2479	q_undi(i,k) = q_seri(i,k)
2480	endif
2481	enddo
2482	enddo
2483
2484	!cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2485	!cc-- pour le soulevement des particules dans le modele convectif
2486	!c
2487	do i = 1,klon
2488	ALE(i) = 0.
2489	ALP(i) = 0.
2490	enddo
2491	!c
2492	!ccalcul de ale_wake et alp_wake
2493	if (iflag_wake>=1) then
2494	if (itap .le. it_wape_prescr) then
2495	do i = 1,klon
2496	ale_wake(i) = wape_prescr
2497	alp_wake(i) = fip_prescr
2498	enddo
2499	else
2500	do i = 1,klon
2501	!cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2502	!ccc ale_wake(i) = 0.5wake_cstar(i)*2
2503	ale_wake(i) = wake_pe(i)
2504	alp_wake(i) = wake_fip(i)
2505	enddo
2506	endif
2507	else
2508	do i = 1,klon
2509	ale_wake(i) = 0.
2510	alp_wake(i) = 0.
2511	enddo
2512	endif
2513	!ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2514	!cdans le thermique sinon
2515	if (iflag_coupl.eq.0) then
2516	if (debut.and.prt_level.gt.9) &
2517	& WRITE(lunout,*)'ALE et ALP imposes'
2518	do i = 1,klon
2519	!con ne couple que ale
2520	!c ALE(i) = max(ale_wake(i),Ale_bl(i))
2521	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2522	!con ne couple que alp
2523	!c ALP(i) = alp_wake(i) + Alp_bl(i)
2524	ALP(i) = alp_wake(i) + alp_bl_prescr
2525	enddo
2526	else
2527	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2528	! do i = 1,klon
2529	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2530	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2531	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2532	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2533	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2534	! enddo
2535
2536	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2537	! Modif FH 2010/04/27. Sans doute temporaire.
2538	! Deux options pour le alp_offset : constant si >?? 0 ou proportionnel ??a
2539	! w si <0
2540	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2541	do i = 1,klon
2542	ALE(i) = max(ale_wake(i),Ale_bl(i))
2543	!ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
2544	if (iflag_trig_bl.ge.1) then
2545	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2546	endif
2547	!ccc fin nrlmd le 10/04/2012
2548	if (alp_offset>=0.) then
2549	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2550	else
2551	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2552	if (alp(i)<0.) then
2553	print*,'ALP ',alp(i),alp_wake(i) &
2554	& ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2555	endif
2556	endif
2557	enddo
2558	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2559
2560	endif
2561	do i=1,klon
2562	if (alp(i)>alp_max) then
2563	IF(prt_level>9)WRITE(lunout,*) &
2564	& 'WARNING SUPER ALP (seuil=',alp_max, &
2565	& '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2566	alp(i)=alp_max
2567	endif
2568	if (ale(i)>ale_max) then
2569	IF(prt_level>9)WRITE(lunout,*) &
2570	& 'WARNING SUPER ALE (seuil=',ale_max, &
2571	& '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2572	ale(i)=ale_max
2573	endif
2574	enddo
2575
2576	!cfin calcul ale et alp
2577	!c=================================================================================================
2578
2579
2580	!c sb, oct02:
2581	!c Schema de convection modularise et vectorise:
2582	!c (driver commun aux versions 3 et 4)
2583	!c
2584	IF (ok_cvl) THEN ! new driver for convectL
2585
2586	IF (type_trac == 'repr') THEN
2587	nbtr_tmp=ntra
2588	ELSE
2589	nbtr_tmp=nbtr
2590	END IF
2591	CALL concvl (iflag_con,iflag_clos, &
2592	& dtime,paprs,pplay,t_undi,q_undi, &
2593	& t_wake,q_wake,wake_s, &
2594	& u_seri,v_seri,tr_seri,nbtr_tmp, &
2595	& ALE,ALP, &
2596	& ema_work1,ema_work2, &
2597	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2598	& rain_con, snow_con, ibas_con, itop_con, sigd, &
2599	& ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2600	& Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2601	& pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2602	! RomP >>>
2603	!! . pmflxr,pmflxs,da,phi,mp,
2604	!! . ftd,fqd,lalim_conv,wght_th) &
2605	& pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2606	& ftd,fqd,lalim_conv,wght_th, &
2607	& ev, ep,epmlmMm,eplaMm, &
2608	& wdtrainA,wdtrainM)
2609	! RomP <<<
2610
2611	!IM begin
2612	!c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2613	!c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2614	!IM end
2615	!IM cf. FH
2616	clwcon0=qcondc
2617	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2618
2619	do i = 1, klon
2620	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2621	enddo
2622
2623	ELSE ! ok_cvl
2624
2625	!c MAF conema3 ne contient pas les traceurs
2626	CALL conema3 (dtime, &
2627	& paprs,pplay,t_seri,q_seri, &
2628	& u_seri,v_seri,tr_seri,ntra, &
2629	& ema_work1,ema_work2, &
2630	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2631	& rain_con, snow_con, ibas_con, itop_con, &
2632	& upwd,dnwd,dnwd0,bas,top, &
2633	& Ma,cape,tvp,rflag, &
2634	& pbase &
2635	& ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2636	& ,clwcon0)
2637
2638	ENDIF ! ok_cvl
2639
2640	!c
2641	!c Correction precip
2642	rain_con = rain_con * cvl_corr
2643	snow_con = snow_con * cvl_corr
2644	!c
2645
2646	IF (.NOT. ok_gust) THEN
2647	do i = 1, klon
2648	wd(i)=0.0
2649	enddo
2650	ENDIF
2651
2652	!c =================================================================== c
2653	!c Calcul des proprietes des nuages convectifs
2654	!c
2655
2656	!c calcul des proprietes des nuages convectifs
2657	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2658	call clouds_gno &
2659	& (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2660
2661	!c =================================================================== c
2662
2663	DO i = 1, klon
2664	itop_con(i) = min(max(itop_con(i),1),klev)
2665	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2666	ENDDO
2667
2668	DO i = 1, klon
2669	ema_pcb(i) = paprs(i,ibas_con(i))
2670	ENDDO
2671	DO i = 1, klon
2672	! L'idicage de itop_con peut cacher un pb potentiel
2673	! FH sous la dictee de JYG, CR
2674	ema_pct(i) = paprs(i,itop_con(i)+1)
2675
2676	if (itop_con(i).gt.klev-3) then
2677	if(prt_level >= 9) then
2678	write(lunout,*)'La convection monte trop haut '
2679	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2680	endif
2681	endif
2682	ENDDO
2683	ELSE IF (iflag_con.eq.0) THEN
2684	write(lunout,*) 'On n appelle pas la convection'
2685	clwcon0=0.
2686	rnebcon0=0.
2687	d_t_con=0.
2688	d_q_con=0.
2689	d_u_con=0.
2690	d_v_con=0.
2691	rain_con=0.
2692	snow_con=0.
2693	bas=1
2694	top=1
2695	ELSE
2696	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2697	CALL abort
2698	ENDIF
2699
2700	!c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2701	!c . d_u_con, d_v_con)
2702
2703	!-----------------------------------------------------------------------------------------
2704	! ajout des tendances de la diffusion turbulente
2705	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2706	!-----------------------------------------------------------------------------------------
2707
2708	if (mydebug) then
2709	call writefield_phy('u_seri',u_seri,llm)
2710	call writefield_phy('v_seri',v_seri,llm)
2711	call writefield_phy('t_seri',t_seri,llm)
2712	call writefield_phy('q_seri',q_seri,llm)
2713	endif
2714
2715	!IM
2716	IF (ip_ebil_phy.ge.2) THEN
2717	ztit='after convect'
2718	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2719	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2720	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2721	call diagphy(airephy,ztit,ip_ebil_phy &
2722	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2723	& , zero_v, rain_con, snow_con, ztsol &
2724	& , d_h_vcol, d_qt, d_ec &
2725	& , fs_bound, fq_bound )
2726	END IF
2727	!C
2728	PRINT *,' Lluis 4 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2729	fname = 'After convection'
2730	varname = 'paprs'
2731	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
2732	varname = 'pplay'
2733	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
2734	varname = 'pphi'
2735	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
2736	varname = 't_seri'
2737	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
2738	varname = 'u_seri'
2739	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
2740	varname = 'v_seri'
2741	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
2742	varname = 'q_seri'
2743	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
2744
2745	IF (check) THEN
2746	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2747	WRITE(lunout,*)"aprescon=", za
2748	zx_t = 0.0
2749	za = 0.0
2750	DO i = 1, klon
2751	za = za + airephy(i)/REAL(klon)
2752	zx_t = zx_t + (rain_con(i)+ &
2753	& snow_con(i))*airephy(i)/REAL(klon)
2754	ENDDO
2755	zx_t = zx_t/za*dtime
2756	WRITE(lunout,*)"Precip=", zx_t
2757	ENDIF
2758	IF (zx_ajustq) THEN
2759	DO i = 1, klon
2760	z_apres(i) = 0.0
2761	ENDDO
2762	DO k = 1, klev
2763	DO i = 1, klon
2764	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2765	& *(paprs(i,k)-paprs(i,k+1))/RG
2766	ENDDO
2767	ENDDO
2768	DO i = 1, klon
2769	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2770	& /z_apres(i)
2771	ENDDO
2772	DO k = 1, klev
2773	DO i = 1, klon
2774	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2775	& z_factor(i).LT.(1.0-1.0E-08)) THEN
2776	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2777	ENDIF
2778	ENDDO
2779	ENDDO
2780	ENDIF
2781	zx_ajustq=.FALSE.
2782
2783	!c
2784	!c=============================================================================
2785	!cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2786	!cpour la couche limite diffuse pour l instant
2787	!c
2788	if (iflag_wake>=1) then
2789	DO k=1,klev
2790	DO i=1,klon
2791	dt_dwn(i,k) = ftd(i,k)
2792	wdt_PBL(i,k) = 0.
2793	dq_dwn(i,k) = fqd(i,k)
2794	wdq_PBL(i,k) = 0.
2795	M_dwn(i,k) = dnwd0(i,k)
2796	M_up(i,k) = upwd(i,k)
2797	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2798	udt_PBL(i,k) = 0.
2799	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2800	udq_PBL(i,k) = 0.
2801	ENDDO
2802	ENDDO
2803
2804	if (iflag_wake==2) then
2805	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2806	DO k = 1,klev
2807	dt_dwn(:,k)= dt_dwn(:,k)+ &
2808	& ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2809	dq_dwn(:,k)= dq_dwn(:,k)+ &
2810	& ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2811	ENDDO
2812	endif
2813	!c
2814	!ccalcul caracteristiques de la poche froide
2815	call calWAKE (paprs,pplay,dtime &
2816	& ,t_seri,q_seri,omega &
2817	& ,dt_dwn,dq_dwn,M_dwn,M_up &
2818	& ,dt_a,dq_a,sigd &
2819	& ,wdt_PBL,wdq_PBL &
2820	& ,udt_PBL,udq_PBL &
2821	& ,wake_deltat,wake_deltaq,wake_dth &
2822	& ,wake_h,wake_s,wake_dens &
2823	& ,wake_pe,wake_fip,wake_gfl &
2824	& ,dt_wake,dq_wake &
2825	& ,wake_k, t_undi,q_undi &
2826	& ,wake_omgbdth,wake_dp_omgb &
2827	& ,wake_dtKE,wake_dqKE &
2828	& ,wake_dtPBL,wake_dqPBL &
2829	& ,wake_omg,wake_dp_deltomg &
2830	& ,wake_spread,wake_Cstar,wake_d_deltat_gw &
2831	& ,wake_ddeltat,wake_ddeltaq)
2832	!c
2833	!-----------------------------------------------------------------------------------------
2834	! ajout des tendances des poches froides
2835	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2836	! coherent avec les autres d_t_...
2837	d_t_wake(:,:)=dt_wake(:,:)*dtime
2838	d_q_wake(:,:)=dq_wake(:,:)*dtime
2839	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2840	!-----------------------------------------------------------------------------------------
2841
2842	endif
2843	!c
2844	!c===================================================================
2845	!cJYG
2846	IF (ip_ebil_phy.ge.2) THEN
2847	ztit='after wake'
2848	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2849	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2850	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2851	call diagphy(airephy,ztit,ip_ebil_phy &
2852	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2853	& , zero_v, zero_v, zero_v, ztsol &
2854	& , d_h_vcol, d_qt, d_ec &
2855	& , fs_bound, fq_bound )
2856	END IF
2857	PRINT *,' Lluis 5 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2858	fname = 'After wakes'
2859	varname = 'paprs'
2860	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
2861	varname = 'pplay'
2862	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
2863	varname = 'pphi'
2864	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
2865	varname = 't_seri'
2866	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
2867	varname = 'u_seri'
2868	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
2869	varname = 'v_seri'
2870	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
2871	varname = 'q_seri'
2872	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
2873
2874	!c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2875	!c
2876	!c===================================================================
2877	!c Convection seche (thermiques ou ajustement)
2878	!c===================================================================
2879	!c
2880	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2881	& ,seuil_inversion,weak_inversion,dthmin)
2882
2883
2884
2885	d_t_ajsb(:,:)=0.
2886	d_q_ajsb(:,:)=0.
2887	d_t_ajs(:,:)=0.
2888	d_u_ajs(:,:)=0.
2889	d_v_ajs(:,:)=0.
2890	d_q_ajs(:,:)=0.
2891	clwcon0th(:,:)=0.
2892	!c
2893	!c fm_therm(:,:)=0.
2894	!c entr_therm(:,:)=0.
2895	!c detr_therm(:,:)=0.
2896	!c
2897
2898	IF(prt_level>9)WRITE(lunout,*) &
2899	& 'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2900	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2901	if(iflag_thermals.lt.0) then
2902	!c Rien
2903	!c ====
2904	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2905
2906
2907	else
2908
2909	!c Thermiques
2910	!c ==========
2911	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2912	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2913
2914
2915	!ccc nrlmd le 10/04/2012
2916	DO k=1,klev+1
2917	DO i=1,klon
2918	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2919	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2920	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2921	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2922	ENDDO
2923	ENDDO
2924	!ccc fin nrlmd le 10/04/2012
2925
2926	if (iflag_thermals>=1) then
2927	call calltherm(pdtphys &
2928	& ,pplay,paprs,pphi,weak_inversion &
2929	& ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
2930	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2931	& ,fm_therm,entr_therm,detr_therm &
2932	& ,zqasc,clwcon0th,lmax_th,ratqscth &
2933	& ,ratqsdiff,zqsatth &
2934	!con rajoute ale et alp, et les caracteristiques de la couche alim &
2935	& ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2936	& ,ztv,zpspsk,ztla,zthl &
2937	!ccc nrlmd le 10/04/2012 &
2938	& ,pbl_tke_input,pctsrf,omega,airephy &
2939	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2940	& ,n2,s2,ale_bl_stat &
2941	& ,therm_tke_max,env_tke_max &
2942	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2943	& ,alp_bl_conv,alp_bl_stat &
2944	!ccc fin nrlmd le 10/04/2012 &
2945	& ,zqla,ztva )
2946
2947	!ccc nrlmd le 10/04/2012
2948	!c-----------Stochastic triggering-----------
2949	if (iflag_trig_bl.ge.1) then
2950	!c
2951	IF (prt_level .GE. 10) THEN
2952	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2953	& cin, ale_bl_stat, alp_bl_stat
2954	ENDIF
2955
2956	!c----Initialisations
2957	do i=1,klon
2958	proba_notrig(i)=1.
2959	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2960	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2961	tau_trig(i)=tau_trig_shallow
2962	else
2963	tau_trig(i)=tau_trig_deep
2964	endif
2965	enddo
2966	!c
2967	IF (prt_level .GE. 10) THEN
2968	print *,'random_notrig, tau_trig ', &
2969	& random_notrig, tau_trig
2970	print *,'s_trig,s2,n2 ', &
2971	& s_trig,s2,n2
2972	ENDIF
2973
2974	!c----Tirage al\'eatoire et calcul de ale_bl_trig
2975	do i=1,klon
2976	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2977	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2978	& (n2(i)*dtime/tau_trig(i))
2979	!c print *, 'proba_notrig(i) ',proba_notrig(i)
2980	if (random_notrig(i) .ge. proba_notrig(i)) then
2981	ale_bl_trig(i)=ale_bl_stat(i)
2982	else
2983	ale_bl_trig(i)=0.
2984	endif
2985	else
2986	proba_notrig(i)=1.
2987	random_notrig(i)=0.
2988	ale_bl_trig(i)=0.
2989	endif
2990	enddo
2991	!c
2992	IF (prt_level .GE. 10) THEN
2993	print *,'proba_notrig, ale_bl_trig ', &
2994	& proba_notrig, ale_bl_trig
2995	ENDIF
2996
2997	endif !(iflag_trig_bl)
2998
2999	!c-----------Statistical closure-----------
3000	if (iflag_clos_bl.ge.1) then
3001
3002	do i=1,klon
3003	alp_bl(i)=alp_bl_stat(i)
3004	enddo
3005
3006	else
3007
3008	alp_bl_stat(:)=0.
3009
3010	endif !(iflag_clos_bl)
3011
3012	IF (prt_level .GE. 10) THEN
3013	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
3014	ENDIF
3015
3016	!ccc fin nrlmd le 10/04/2012
3017
3018	! ----------------------------------------------------------------------
3019	! Transport de la TKE par les panaches thermiques.
3020	! FH : 2010/02/01
3021	! if (iflag_pbl.eq.10) then
3022	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3023	! s rg,paprs,pbl_tke)
3024	! endif
3025	! ----------------------------------------------------------------------
3026	!IM/FH: 2011/02/23
3027	! Couplage Thermiques/Emanuel seulement si T<0
3028	if (iflag_coupl==2) then
3029	print*,'Couplage Thermiques/Emanuel seulement si T<0'
3030	do i=1,klon
3031	if (t_seri(i,lmax_th(i))>273.) then
3032	Ale_bl(i)=0.
3033	endif
3034	enddo
3035	endif
3036
3037	do i=1,klon
3038	zmax_th(i)=pphi(i,lmax_th(i))/rg
3039	enddo
3040
3041	endif
3042
3043
3044	!c Ajustement sec
3045	!c ==============
3046
3047	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3048	! a partir du sommet des thermiques.
3049	! Dans le cas contraire, on demarre au niveau 1.
3050
3051	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
3052
3053	if(iflag_thermals.eq.0) then
3054	IF(prt_level>9)WRITE(lunout,*)'ajsec'
3055	limbas(:)=1
3056	else
3057	limbas(:)=lmax_th(:)
3058	endif
3059
3060	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3061	! pour des test de convergence numerique.
3062	! Le nouveau ajsec est a priori mieux, meme pour le cas
3063	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3064	! non nulles numeriquement pour des mailles non concernees.
3065
3066	if (iflag_thermals.eq.0) then
3067	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3068	& , d_t_ajsb, d_q_ajsb)
3069	else
3070	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3071	& , d_t_ajsb, d_q_ajsb)
3072	endif
3073
3074	!-----------------------------------------------------------------------------------------
3075	! ajout des tendances de l'ajustement sec ou des thermiques
3076	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
3077	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3078	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3079
3080	!-----------------------------------------------------------------------------------------
3081
3082	endif
3083
3084	endif
3085
3086	!c
3087	!c===================================================================
3088	!IM
3089	IF (ip_ebil_phy.ge.2) THEN
3090	ztit='after dry_adjust'
3091	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3092	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3093	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3094	call diagphy(airephy,ztit,ip_ebil_phy &
3095	& , zero_v, zero_v, zero_v, zero_v, zero_v &
3096	& , zero_v, zero_v, zero_v, ztsol &
3097	& , d_h_vcol, d_qt, d_ec &
3098	& , fs_bound, fq_bound )
3099	END IF
3100	PRINT *,' Lluis 6 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3101	fname = 'After thermals'
3102	varname = 'paprs'
3103	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
3104	varname = 'pplay'
3105	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
3106	varname = 'pphi'
3107	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
3108	varname = 't_seri'
3109	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
3110	varname = 'u_seri'
3111	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
3112	varname = 'v_seri'
3113	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
3114	varname = 'q_seri'
3115	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
3116
3117	!c-------------------------------------------------------------------------
3118	! Computation of ratqs, the width (normalized) of the subrid scale
3119	! water distribution
3120	CALL calcratqs(klon,klev,prt_level,lunout, &
3121	& iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
3122	& ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
3123	& ptconv,ptconvth,clwcon0th, rnebcon0th, &
3124	& paprs,pplay,q_seri,zqsat,fm_therm, &
3125	& ratqs,ratqsc)
3126
3127
3128	!c
3129	!c Appeler le processus de condensation a grande echelle
3130	!c et le processus de precipitation
3131	!c-------------------------------------------------------------------------
3132	IF (prt_level .GE.10) THEN
3133	print *,' ->fisrtilp '
3134	ENDIF
3135	!c-------------------------------------------------------------------------
3136	CALL fisrtilp(dtime,paprs,pplay, &
3137	& t_seri, q_seri,ptconv,ratqs, &
3138	& d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
3139	& rain_lsc, snow_lsc, &
3140	& pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3141	& frac_impa, frac_nucl, beta_prec_fisrt, &
3142	& prfl, psfl, rhcl, &
3143	& zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
3144
3145	WHERE (rain_lsc < 0) rain_lsc = 0.
3146	WHERE (snow_lsc < 0) snow_lsc = 0.
3147	!-----------------------------------------------------------------------------------------
3148	! ajout des tendances de la diffusion turbulente
3149	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
3150	!-----------------------------------------------------------------------------------------
3151	DO k = 1, klev
3152	DO i = 1, klon
3153	cldfra(i,k) = rneb(i,k)
3154	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3155	ENDDO
3156	ENDDO
3157	IF (check) THEN
3158	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
3159	WRITE(lunout,*)"apresilp=", za
3160	zx_t = 0.0
3161	za = 0.0
3162	DO i = 1, klon
3163	za = za + airephy(i)/REAL(klon)
3164	zx_t = zx_t + (rain_lsc(i) &
3165	& + snow_lsc(i))*airephy(i)/REAL(klon)
3166	ENDDO
3167	zx_t = zx_t/za*dtime
3168	WRITE(lunout,*)"Precip=", zx_t
3169	ENDIF
3170	!IM
3171	IF (ip_ebil_phy.ge.2) THEN
3172	ztit='after fisrt'
3173	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3174	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3175	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3176	call diagphy(airephy,ztit,ip_ebil_phy &
3177	& , zero_v, zero_v, zero_v, zero_v, zero_v &
3178	& , zero_v, rain_lsc, snow_lsc, ztsol &
3179	& , d_h_vcol, d_qt, d_ec &
3180	& , fs_bound, fq_bound )
3181	END IF
3182	PRINT *,' Lluis 7 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3183	fname = 'After firstilp'
3184	varname = 'paprs'
3185	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
3186	varname = 'pplay'
3187	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
3188	varname = 'pphi'
3189	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
3190	varname = 't_seri'
3191	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
3192	varname = 'u_seri'
3193	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
3194	varname = 'v_seri'
3195	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
3196	varname = 'q_seri'
3197	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
3198
3199	if (mydebug) then
3200	call writefield_phy('u_seri',u_seri,llm)
3201	call writefield_phy('v_seri',v_seri,llm)
3202	call writefield_phy('t_seri',t_seri,llm)
3203	call writefield_phy('q_seri',q_seri,llm)
3204	endif
3205	!c
3206	!c-------------------------------------------------------------------
3207	!c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3208	!c-------------------------------------------------------------------
3209
3210	!c 1. NUAGES CONVECTIFS
3211	!c
3212	!IM cf FH
3213	!c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
3214	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
3215	snow_tiedtke=0.
3216	!c print*,'avant calcul de la pseudo precip '
3217	!c print*,'iflag_cldcon',iflag_cldcon
3218	if (iflag_cldcon.eq.-1) then
3219	rain_tiedtke=rain_con
3220	else
3221	!c print*,'calcul de la pseudo precip '
3222	rain_tiedtke=0.
3223	!c print*,'calcul de la pseudo precip 0'
3224	do k=1,klev
3225	do i=1,klon
3226	if (d_q_con(i,k).lt.0.) then
3227	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3228	& *(paprs(i,k)-paprs(i,k+1))/rg
3229	endif
3230	enddo
3231	enddo
3232	endif
3233	!c
3234	!c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3235	!c
3236
3237	!c Nuages diagnostiques pour Tiedtke
3238	CALL diagcld1(paprs,pplay, &
3239	!IM cf FH . rain_con,snow_con,ibas_con,itop_con, &
3240	& rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3241	& diafra,dialiq)
3242	DO k = 1, klev
3243	DO i = 1, klon
3244	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3245	cldliq(i,k) = dialiq(i,k)
3246	cldfra(i,k) = diafra(i,k)
3247	ENDIF
3248	ENDDO
3249	ENDDO
3250
3251	ELSE IF (iflag_cldcon.ge.3) THEN
3252	!c On prend pour les nuages convectifs le max du calcul de la
3253	!c convection et du calcul du pas de temps precedent diminue d'un facteur
3254	!c facttemps
3255	facteur = pdtphys *facttemps
3256	do k=1,klev
3257	do i=1,klon
3258	rnebcon(i,k)=rnebcon(i,k)*facteur
3259	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
3260	& then
3261	rnebcon(i,k)=rnebcon0(i,k)
3262	clwcon(i,k)=clwcon0(i,k)
3263	endif
3264	enddo
3265	enddo
3266
3267	!c
3268	!cjq - introduce the aerosol direct and first indirect radiative forcings
3269	!cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3270	IF (flag_aerosol .gt. 0) THEN
3271	IF (.NOT. aerosol_couple) &
3272	& CALL readaerosol_optic( &
3273	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3274	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3275	& mass_solu_aero, mass_solu_aero_pi, &
3276	& tau_aero, piz_aero, cg_aero, &
3277	& tausum_aero, tau3d_aero)
3278	ELSE
3279	tausum_aero(:,:,:) = 0.
3280	tau_aero(:,:,:,:) = 0.
3281	piz_aero(:,:,:,:) = 0.
3282	cg_aero(:,:,:,:) = 0.
3283	ENDIF
3284	!c
3285	!c--STRAT AEROSOL
3286	!c--updates tausum_aero,tau_aero,piz_aero,cg_aero
3287	IF (flag_aerosol_strat) THEN
3288	PRINT *,'appel a readaerosolstrat', mth_cur
3289	CALL readaerosolstrato(debut)
3290	ENDIF
3291	!c--fin STRAT AEROSOL
3292
3293	!IM calcul nuages par le simulateur ISCCP
3294	!c
3295	#ifdef histISCCP
3296	IF (ok_isccp) THEN
3297	!c
3298	!IM lecture invtau, tautab des fichiers formattes
3299	!c
3300	IF (debut) THEN
3301	!$OMP MASTER
3302	!c
3303	open(99,file='tautab.formatted', FORM='FORMATTED')
3304	read(99,'(f30.20)') tautab_omp
3305	close(99)
3306	!c
3307	open(99,file='invtau.formatted',form='FORMATTED')
3308	read(99,'(i10)') invtau_omp
3309
3310	!c print*,'calcul_simulISCCP invtau_omp',invtau_omp
3311	!c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
3312
3313	close(99)
3314	!$OMP END MASTER
3315	!$OMP BARRIER
3316	tautab=tautab_omp
3317	invtau=invtau_omp
3318	!c
3319	ENDIF !debut
3320	!c
3321	!IM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
3322	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
3323	#include "calcul_simulISCCP.h"
3324	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
3325	ENDIF !ok_isccp
3326	#endif
3327
3328	!c On prend la somme des fractions nuageuses et des contenus en eau
3329
3330	if (iflag_cldcon>=5) then
3331
3332	do k=1,klev
3333	ptconvth(:,k)=fm_therm(:,k+1)>0.
3334	enddo
3335
3336	if (iflag_coupl==4) then
3337
3338	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3339	! convectives et lsc dans la partie des thermiques
3340	! Le controle par iflag_coupl est peut etre provisoire.
3341	do k=1,klev
3342	do i=1,klon
3343	if (ptconv(i,k).and.ptconvth(i,k)) then
3344	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3345	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3346	else if (ptconv(i,k)) then
3347	cldfra(i,k)=rnebcon(i,k)
3348	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3349	endif
3350	enddo
3351	enddo
3352
3353	else if (iflag_coupl==5) then
3354	do k=1,klev
3355	do i=1,klon
3356	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3357	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3358	enddo
3359	enddo
3360
3361	else
3362
3363	! Si on est sur un point touche par la convection profonde et pas
3364	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3365	! de la convection profonde.
3366
3367	!IM/FH: 2011/02/23
3368	! definition des points sur lesquels ls thermiques sont actifs
3369
3370	do k=1,klev
3371	do i=1,klon
3372	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3373	cldfra(i,k)=rnebcon(i,k)
3374	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3375	endif
3376	enddo
3377	enddo
3378
3379	endif
3380
3381	else
3382
3383	! Ancienne version
3384	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3385	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3386	endif
3387
3388	ENDIF
3389
3390	! plulsc(:)=0.
3391	! do k=1,klev,-1
3392	! do i=1,klon
3393	! zzz=prfl(:,k)+psfl(:,k)
3394	! if (.not.ptconvth.zzz.gt.0.)
3395	! enddo prfl, psfl,
3396	! enddo
3397	!c
3398	!c 2. NUAGES STARTIFORMES
3399	!c
3400	IF (ok_stratus) THEN
3401	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3402	DO k = 1, klev
3403	DO i = 1, klon
3404	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3405	cldliq(i,k) = dialiq(i,k)
3406	cldfra(i,k) = diafra(i,k)
3407	ENDIF
3408	ENDDO
3409	ENDDO
3410	ENDIF
3411	!c
3412	!c Precipitation totale
3413	!c
3414	DO i = 1, klon
3415	rain_fall(i) = rain_con(i) + rain_lsc(i)
3416	snow_fall(i) = snow_con(i) + snow_lsc(i)
3417	ENDDO
3418	!IM
3419	IF (ip_ebil_phy.ge.2) THEN
3420	ztit="after diagcld"
3421	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3422	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3423	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3424	call diagphy(airephy,ztit,ip_ebil_phy &
3425	& , zero_v, zero_v, zero_v, zero_v, zero_v &
3426	& , zero_v, zero_v, zero_v, ztsol &
3427	& , d_h_vcol, d_qt, d_ec &
3428	& , fs_bound, fq_bound )
3429	END IF
3430	PRINT *,' Lluis 8 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3431	fname = 'After coupling convection+thermics+wakes'
3432	varname = 'paprs'
3433	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
3434	varname = 'pplay'
3435	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
3436	varname = 'pphi'
3437	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
3438	varname = 't_seri'
3439	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
3440	varname = 'u_seri'
3441	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
3442	varname = 'v_seri'
3443	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
3444	varname = 'q_seri'
3445	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
3446
3447	!c
3448	!c Calculer l'humidite relative pour diagnostique
3449	!c
3450	DO k = 1, klev
3451	DO i = 1, klon
3452	zx_t = t_seri(i,k)
3453	IF (thermcep) THEN
3454	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3455	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3456	zx_qs = MIN(0.5,zx_qs)
3457	zcor = 1./(1.-retv*zx_qs)
3458	zx_qs = zx_qs*zcor
3459	ELSE
3460	IF (zx_t.LT.t_coup) THEN
3461	zx_qs = qsats(zx_t)/pplay(i,k)
3462	ELSE
3463	zx_qs = qsatl(zx_t)/pplay(i,k)
3464	ENDIF
3465	ENDIF
3466	zx_rh(i,k) = q_seri(i,k)/zx_qs
3467	zqsat(i,k)=zx_qs
3468	ENDDO
3469	ENDDO
3470
3471	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3472	!c equivalente a 2m (tpote) pour diagnostique
3473	!c
3474	DO i = 1, klon
3475	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3476	IF (thermcep) THEN
3477	IF(zt2m(i).LT.RTT) then
3478	Lheat=RLSTT
3479	ELSE
3480	Lheat=RLVTT
3481	ENDIF
3482	ELSE
3483	IF (zt2m(i).LT.RTT) THEN
3484	Lheat=RLSTT
3485	ELSE
3486	Lheat=RLVTT
3487	ENDIF
3488	ENDIF
3489	tpote(i) = tpot(i)* &
3490	& EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3491	ENDDO
3492
3493	IF (type_trac == 'inca') THEN
3494	#ifdef INCA
3495	CALL VTe(VTphysiq)
3496	CALL VTb(VTinca)
3497	calday = REAL(days_elapsed + 1) + jH_cur
3498
3499	call chemtime(itap+itau_phy-1, date0, dtime)
3500	IF (config_inca == 'aero') THEN
3501	CALL AEROSOL_METEO_CALC( &
3502	& calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3503	& prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3504	END IF
3505
3506	zxsnow_dummy(:) = 0.0
3507
3508	CALL chemhook_begin (calday, &
3509	& days_elapsed+1, &
3510	& jH_cur, &
3511	& pctsrf(1,1), &
3512	& rlat, &
3513	& rlon, &
3514	& airephy, &
3515	& paprs, &
3516	& pplay, &
3517	& coefh(:,:,is_ave), &
3518	& pphi, &
3519	& t_seri, &
3520	& u, &
3521	& v, &
3522	& wo(:, :, 1), &
3523	& q_seri, &
3524	& zxtsol, &
3525	& zxsnow_dummy, &
3526	& solsw, &
3527	& albsol1, &
3528	& rain_fall, &
3529	& snow_fall, &
3530	& itop_con, &
3531	& ibas_con, &
3532	& cldfra, &
3533	& iim, &
3534	& jjm, &
3535	& tr_seri, &
3536	& ftsol, &
3537	& paprs, &
3538	& cdragh, &
3539	& cdragm, &
3540	& pctsrf, &
3541	& pdtphys, &
3542	& itap)
3543
3544	CALL VTe(VTinca)
3545	CALL VTb(VTphysiq)
3546	#endif
3547	END IF !type_trac = inca
3548	!c
3549	!c Calculer les parametres optiques des nuages et quelques
3550	!c parametres pour diagnostiques:
3551	!c
3552
3553	IF (aerosol_couple) THEN
3554	mass_solu_aero(:,:) = ccm(:,:,1)
3555	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3556	END IF
3557
3558	if (ok_newmicro) then
3559	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3560	& paprs, pplay, t_seri, cldliq, cldfra, &
3561	& cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3562	& flwp, fiwp, flwc, fiwc, &
3563	& mass_solu_aero, mass_solu_aero_pi, &
3564	& cldtaupi, re, fl, ref_liq, ref_ice)
3565	else
3566	CALL nuage (paprs, pplay, &
3567	& t_seri, cldliq, cldfra, cldtau, cldemi, &
3568	& cldh, cldl, cldm, cldt, cldq, &
3569	& ok_aie, &
3570	& mass_solu_aero, mass_solu_aero_pi, &
3571	& bl95_b0, bl95_b1, &
3572	& cldtaupi, re, fl)
3573	endif
3574	!c
3575	!IM betaCRF
3576	!c
3577	cldtaurad = cldtau
3578	cldtaupirad = cldtaupi
3579	cldemirad = cldemi &
3580	&!c
3581	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3582	&lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3583	!c
3584	!c global
3585	!c
3586	DO k=1, klev
3587	DO i=1, klon
3588	if (pplay(i,k).GE.pfree) THEN
3589	beta(i,k) = beta_pbl
3590	else
3591	beta(i,k) = beta_free
3592	endif
3593	if (mskocean_beta) THEN
3594	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3595	endif
3596	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3597	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3598	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3599	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3600	ENDDO
3601	ENDDO
3602	!c
3603	else
3604	!c
3605	!c regional
3606	!c
3607	DO k=1, klev
3608	DO i=1,klon
3609	!c
3610	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3611	& rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3612	if (pplay(i,k).GE.pfree) THEN
3613	beta(i,k) = beta_pbl
3614	else
3615	beta(i,k) = beta_free
3616	endif
3617	if (mskocean_beta) THEN
3618	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3619	endif
3620	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3621	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3622	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3623	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3624	endif
3625	!c
3626	ENDDO
3627	ENDDO
3628	!c
3629	endif
3630	!c
3631	!c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3632	!c
3633	IF (MOD(itaprad,radpas).EQ.0) THEN
3634
3635	DO i = 1, klon
3636	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3637	& + falb1(i,is_lic) * pctsrf(i,is_lic) &
3638	& + falb1(i,is_ter) * pctsrf(i,is_ter) &
3639	& + falb1(i,is_sic) * pctsrf(i,is_sic)
3640	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3641	& + falb2(i,is_lic) * pctsrf(i,is_lic) &
3642	& + falb2(i,is_ter) * pctsrf(i,is_ter) &
3643	& + falb2(i,is_sic) * pctsrf(i,is_sic)
3644	ENDDO
3645
3646	if (mydebug) then
3647	call writefield_phy('u_seri',u_seri,llm)
3648	call writefield_phy('v_seri',v_seri,llm)
3649	call writefield_phy('t_seri',t_seri,llm)
3650	call writefield_phy('q_seri',q_seri,llm)
3651	endif
3652
3653	IF (aerosol_couple) THEN
3654	#ifdef INCA
3655	CALL radlwsw_inca &
3656	& (kdlon,kflev,dist, rmu0, fract, solaire, &
3657	& paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3658	& wo(:, :, 1), &
3659	& cldfrarad, cldemirad, cldtaurad, &
3660	& heat,heat0,cool,cool0,radsol,albpla, &
3661	& topsw,toplw,solsw,sollw, &
3662	& sollwdown, &
3663	& topsw0,toplw0,solsw0,sollw0, &
3664	& lwdn0, lwdn, lwup0, lwup, &
3665	& swdn0, swdn, swup0, swup, &
3666	& ok_ade, ok_aie, &
3667	& tau_aero, piz_aero, cg_aero, &
3668	& topswad_aero, solswad_aero, &
3669	& topswad0_aero, solswad0_aero, &
3670	& topsw_aero, topsw0_aero, &
3671	& solsw_aero, solsw0_aero, &
3672	& cldtaupirad, &
3673	& topswai_aero, solswai_aero)
3674
3675	#endif
3676	ELSE
3677	!c
3678	!IM calcul radiatif pour le cas actuel
3679	!c
3680	RCO2 = RCO2_act
3681	RCH4 = RCH4_act
3682	RN2O = RN2O_act
3683	RCFC11 = RCFC11_act
3684	RCFC12 = RCFC12_act
3685	!c
3686	IF (prt_level .GE.10) THEN
3687	print *,' ->radlwsw, number 1 '
3688	ENDIF
3689	!c
3690	CALL radlwsw &
3691	& (dist, rmu0, fract, &
3692	& paprs, pplay,zxtsol,albsol1, albsol2, &
3693	& t_seri,q_seri,wo, &
3694	& cldfrarad, cldemirad, cldtaurad, &
3695	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3696	& flag_aerosol_strat, &
3697	& tau_aero, piz_aero, cg_aero, &
3698	& cldtaupirad,new_aod, &
3699	& zqsat, flwc, fiwc, &
3700	& heat,heat0,cool,cool0,radsol,albpla, &
3701	& topsw,toplw,solsw,sollw, &
3702	& sollwdown, &
3703	& topsw0,toplw0,solsw0,sollw0, &
3704	& lwdn0, lwdn, lwup0, lwup, &
3705	& swdn0, swdn, swup0, swup, &
3706	& topswad_aero, solswad_aero, &
3707	& topswai_aero, solswai_aero, &
3708	& topswad0_aero, solswad0_aero, &
3709	& topsw_aero, topsw0_aero, &
3710	& solsw_aero, solsw0_aero, &
3711	& topswcf_aero, solswcf_aero)
3712
3713	!c
3714	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3715	!IM des taux doit etre different du taux actuel
3716	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3717	!c
3718
3719	if (ok_4xCO2atm) then
3720	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3721	&RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3722	&RCFC12_per.NE.RCFC12_act) THEN
3723	!c
3724	RCO2 = RCO2_per
3725	RCH4 = RCH4_per
3726	RN2O = RN2O_per
3727	RCFC11 = RCFC11_per
3728	RCFC12 = RCFC12_per
3729	!c
3730	IF (prt_level .GE.10) THEN
3731	print *,' ->radlwsw, number 2 '
3732	ENDIF
3733	!c
3734	CALL radlwsw &
3735	& (dist, rmu0, fract, &
3736	& paprs, pplay,zxtsol,albsol1, albsol2, &
3737	& t_seri,q_seri,wo, &
3738	& cldfra, cldemi, cldtau, &
3739	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3740	& flag_aerosol_strat, &
3741	& tau_aero, piz_aero, cg_aero, &
3742	& cldtaupi,new_aod, &
3743	& zqsat, flwc, fiwc, &
3744	& heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3745	& topswp,toplwp,solswp,sollwp, &
3746	& sollwdownp, &
3747	& topsw0p,toplw0p,solsw0p,sollw0p, &
3748	& lwdn0p, lwdnp, lwup0p, lwupp, &
3749	& swdn0p, swdnp, swup0p, swupp, &
3750	& topswad_aerop, solswad_aerop, &
3751	& topswai_aerop, solswai_aerop, &
3752	& topswad0_aerop, solswad0_aerop, &
3753	& topsw_aerop, topsw0_aerop, &
3754	& solsw_aerop, solsw0_aerop, &
3755	& topswcf_aerop, solswcf_aerop)
3756	endif
3757	endif
3758	!c
3759	ENDIF ! aerosol_couple
3760	itaprad = 0
3761	ENDIF ! MOD(itaprad,radpas)
3762	itaprad = itaprad + 1
3763
3764	IF (iflag_radia.eq.0) THEN
3765	IF (prt_level.ge.9) THEN
3766	PRINT *,'--------------------------------------------------'
3767	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3768	PRINT *,'>>>> heat et cool mis a zero '
3769	PRINT *,'--------------------------------------------------'
3770	END IF
3771	heat=0.
3772	cool=0.
3773	sollw=0. ! MPL 01032011
3774	solsw=0.
3775	radsol=0.
3776	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3777	swup0=0.
3778	swdn=0.
3779	swdn0=0.
3780	lwup=0.
3781	lwup0=0.
3782	lwdn=0.
3783	lwdn0=0.
3784	END IF
3785
3786	!c
3787	!c Ajouter la tendance des rayonnements (tous les pas)
3788	!c
3789	DO k = 1, klev
3790	DO i = 1, klon
3791	t_seri(i,k) = t_seri(i,k) &
3792	& + (heat(i,k)-cool(i,k)) * dtime/RDAY
3793	ENDDO
3794	ENDDO
3795	!c
3796	if (mydebug) then
3797	call writefield_phy('u_seri',u_seri,llm)
3798	call writefield_phy('v_seri',v_seri,llm)
3799	call writefield_phy('t_seri',t_seri,llm)
3800	call writefield_phy('q_seri',q_seri,llm)
3801	endif
3802
3803	!IM
3804	IF (ip_ebil_phy.ge.2) THEN
3805	ztit='after rad'
3806	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3807	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3808	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3809	call diagphy(airephy,ztit,ip_ebil_phy &
3810	& , topsw, toplw, solsw, sollw, zero_v &
3811	& , zero_v, zero_v, zero_v, ztsol &
3812	& , d_h_vcol, d_qt, d_ec &
3813	& , fs_bound, fq_bound )
3814	END IF
3815	PRINT *,' Lluis 9 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3816	fname = 'After radiation'
3817	varname = 'paprs'
3818	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
3819	varname = 'pplay'
3820	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
3821	varname = 'pphi'
3822	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
3823	varname = 't_seri'
3824	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
3825	varname = 'u_seri'
3826	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
3827	varname = 'v_seri'
3828	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
3829	varname = 'q_seri'
3830	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
3831
3832	!c
3833	!c
3834	!c Calculer l'hydrologie de la surface
3835	!c
3836	!c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3837	!c . agesno, ftsol,fqsurf,fsnow, ruis)
3838	!c
3839
3840	!c
3841	!c Calculer le bilan du sol et la derive de temperature (couplage)
3842	!c
3843	DO i = 1, klon
3844	!c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3845	!c a la demande de JLD
3846	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3847	ENDDO
3848	!c
3849	!cmoddeblott(jan95)
3850	!c Appeler le programme de parametrisation de l'orographie
3851	!c a l'echelle sous-maille:
3852	!c
3853	IF (prt_level .GE.10) THEN
3854	print *,' call orography ? ', ok_orodr
3855	ENDIF
3856	!c
3857	IF (ok_orodr) THEN
3858	!c
3859	!c selection des points pour lesquels le shema est actif:
3860	igwd=0
3861	DO i=1,klon
3862	itest(i)=0
3863	!c IF ((zstd(i).gt.10.0)) THEN
3864	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3865	itest(i)=1
3866	igwd=igwd+1
3867	idx(igwd)=i
3868	ENDIF
3869	ENDDO
3870	!c igwdim=MAX(1,igwd)
3871	!c
3872	IF (ok_strato) THEN
3873
3874	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3875	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3876	& igwd,idx,itest, &
3877	& t_seri, u_seri, v_seri, &
3878	& zulow, zvlow, zustrdr, zvstrdr, &
3879	& d_t_oro, d_u_oro, d_v_oro)
3880
3881	ELSE
3882	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3883	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3884	& igwd,idx,itest, &
3885	& t_seri, u_seri, v_seri, &
3886	& zulow, zvlow, zustrdr, zvstrdr, &
3887	& d_t_oro, d_u_oro, d_v_oro)
3888	ENDIF
3889	!c
3890	!c ajout des tendances
3891	!-----------------------------------------------------------------------------------------
3892	! ajout des tendances de la trainee de l'orographie
3893	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3894	!-----------------------------------------------------------------------------------------
3895	!c
3896	ENDIF ! fin de test sur ok_orodr
3897	!c
3898	if (mydebug) then
3899	call writefield_phy('u_seri',u_seri,llm)
3900	call writefield_phy('v_seri',v_seri,llm)
3901	call writefield_phy('t_seri',t_seri,llm)
3902	call writefield_phy('q_seri',q_seri,llm)
3903	endif
3904
3905	IF (ok_orolf) THEN
3906	!c
3907	!c selection des points pour lesquels le shema est actif:
3908	igwd=0
3909	DO i=1,klon
3910	itest(i)=0
3911	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3912	itest(i)=1
3913	igwd=igwd+1
3914	idx(igwd)=i
3915	ENDIF
3916	ENDDO
3917	!c igwdim=MAX(1,igwd)
3918	!c
3919	IF (ok_strato) THEN
3920
3921	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3922	& rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3923	& igwd,idx,itest, &
3924	& t_seri, u_seri, v_seri, &
3925	& zulow, zvlow, zustrli, zvstrli, &
3926	& d_t_lif, d_u_lif, d_v_lif )
3927
3928	ELSE
3929	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3930	& rlat,zmea,zstd,zpic, &
3931	& itest, &
3932	& t_seri, u_seri, v_seri, &
3933	& zulow, zvlow, zustrli, zvstrli, &
3934	& d_t_lif, d_u_lif, d_v_lif)
3935	ENDIF
3936	!c
3937	!-----------------------------------------------------------------------------------------
3938	! ajout des tendances de la portance de l'orographie
3939	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3940	!-----------------------------------------------------------------------------------------
3941	!c
3942	ENDIF ! fin de test sur ok_orolf
3943	!C HINES GWD PARAMETRIZATION
3944
3945	IF (ok_hines) then
3946
3947	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3948	& rlat,t_seri,u_seri,v_seri, &
3949	& zustrhi,zvstrhi, &
3950	& d_t_hin, d_u_hin, d_v_hin)
3951	!c
3952	!c ajout des tendances
3953	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
3954
3955	ENDIF
3956	!c
3957
3958	!c
3959	!IM cf. FLott BEG
3960	!C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3961
3962	if (mydebug) then
3963	call writefield_phy('u_seri',u_seri,llm)
3964	call writefield_phy('v_seri',v_seri,llm)
3965	call writefield_phy('t_seri',t_seri,llm)
3966	call writefield_phy('q_seri',q_seri,llm)
3967	endif
3968
3969	DO i = 1, klon
3970	zustrph(i)=0.
3971	zvstrph(i)=0.
3972	ENDDO
3973	DO k = 1, klev
3974	DO i = 1, klon
3975	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3976	& (paprs(i,k)-paprs(i,k+1))/rg
3977	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3978	& (paprs(i,k)-paprs(i,k+1))/rg
3979	ENDDO
3980	ENDDO
3981	!c
3982	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3983	!c
3984	! IF (is_sequential .and. ok_orodr) THEN
3985	! CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3986	! & ra,rg,romega, &
3987	! & rlat,rlon,pphis, &
3988	! & zustrdr,zustrli,zustrph, &
3989	! & zvstrdr,zvstrli,zvstrph, &
3990	! & paprs,u,v, &
3991	! & aam, torsfc)
3992	! ENDIF
3993	!IM cf. FLott END
3994	!IM
3995	IF (ip_ebil_phy.ge.2) THEN
3996	ztit='after orography'
3997	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3998	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3999	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
4000	call diagphy(airephy,ztit,ip_ebil_phy &
4001	& , zero_v, zero_v, zero_v, zero_v, zero_v &
4002	& , zero_v, zero_v, zero_v, ztsol &
4003	& , d_h_vcol, d_qt, d_ec &
4004	& , fs_bound, fq_bound )
4005	END IF
4006	PRINT *,' Lluis 10 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
4007	fname = 'after orography'
4008	varname = 'paprs'
4009	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
4010	varname = 'pplay'
4011	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
4012	varname = 'pphi'
4013	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
4014	varname = 't_seri'
4015	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
4016	varname = 'u_seri'
4017	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
4018	varname = 'v_seri'
4019	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
4020	varname = 'q_seri'
4021	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
4022	!c
4023	!c
4024	!====================================================================
4025	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4026	!====================================================================
4027	! Abderrahmane 24.08.09
4028
4029	IF (ok_cosp) THEN
4030	! adeclarer
4031	#ifdef CPP_COSP
4032	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4033
4034	print*,'freq_cosp',freq_cosp
4035	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4036	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4037	! s ref_liq,ref_ice
4038	call phys_cosp(itap,dtime,freq_cosp, &
4039	& ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4040	& ecrit_mth,ecrit_day,ecrit_hf, &
4041	& klon,klev,rlon,rlat,presnivs,overlap, &
4042	& ref_liq,ref_ice, &
4043	& pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4044	& zu10m,zv10m,pphis, &
4045	& zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4046	& qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4047	& prfl(:,1:klev),psfl(:,1:klev), &
4048	& pmflxr(:,1:klev),pmflxs(:,1:klev), &
4049	& mr_ozone,cldtau, cldemi)
4050
4051	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4052	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4053	! M clMISR,
4054	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4055	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4056
4057	ENDIF
4058
4059	#endif
4060	ENDIF !ok_cosp
4061	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4062	!cAA
4063	!cAA Installation de l'interface online-offline pour traceurs
4064	!cAA
4065	!c====================================================================
4066	!c Calcul des tendances traceurs
4067	!c====================================================================
4068	!C
4069
4070	IF (type_trac=='repr') THEN
4071	sh_in(:,:) = q_seri(:,:)
4072	ELSE
4073	sh_in(:,:) = qx(:,:,ivap)
4074	END IF
4075
4076	call phytrac ( &
4077	& itap, days_elapsed+1, jH_cur, debut, &
4078	& lafin, dtime, u, v, t, &
4079	& paprs, pplay, pmfu, pmfd, &
4080	& pen_u, pde_u, pen_d, pde_d, &
4081	& cdragh, coefh(:,:,is_ave), fm_therm, entr_therm, &
4082	& u1, v1, ftsol, pctsrf, &
4083	& ustar, u10m, v10m, &
4084	& rlat, rlon, &
4085	& frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4086	& presnivs, pphis, pphi, albsol1, &
4087	& sh_in, rhcl, cldfra, rneb, &
4088	& diafra, cldliq, itop_con, ibas_con, &
4089	& pmflxr, pmflxs, prfl, psfl, &
4090	& da, phi, mp, upwd, &
4091	& phi2, d1a, dam, sij, &
4092	& wdtrainA, wdtrainM, sigd, clw,elij, &
4093	& ev, ep, epmlmMm, eplaMm, &
4094	& dnwd, aerosol_couple, flxmass_w, &
4095	& tau_aero, piz_aero, cg_aero, ccm, &
4096	& rfname, &
4097	& d_tr_dyn, &
4098	& tr_seri)
4099
4100	IF (offline) THEN
4101
4102	IF (prt_level.ge.9) &
4103	& print*,'Attention on met a 0 les thermiques pour phystoke'
4104	call phystokenc ( &
4105	& nlon,klev,pdtphys,rlon,rlat, &
4106	& t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4107	& fm_therm,entr_therm, &
4108	& cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf, &
4109	& frac_impa, frac_nucl, &
4110	& pphis,airephy,dtime,itap, &
4111	& qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4112
4113
4114	ENDIF
4115
4116	!c
4117	!c Calculer le transport de l'eau et de l'energie (diagnostique)
4118	!c
4119	CALL transp (paprs,zxtsol, &
4120	& t_seri, q_seri, u_seri, v_seri, zphi, &
4121	& ve, vq, ue, uq)
4122	!c
4123	!IM global posePB BEG
4124	IF(1.EQ.0) THEN
4125	!c
4126	CALL transp_lay (paprs,zxtsol, &
4127	& t_seri, q_seri, u_seri, v_seri, zphi, &
4128	& ve_lay, vq_lay, ue_lay, uq_lay)
4129	!c
4130	ENDIF !(1.EQ.0) THEN
4131	!IM global posePB END
4132	!c Accumuler les variables a stocker dans les fichiers histoire:
4133	!c
4134	fname = 'after phytrac'
4135	varname = 'paprs'
4136	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
4137	varname = 'pplay'
4138	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
4139	varname = 'pphi'
4140	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
4141	varname = 't_seri'
4142	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
4143	varname = 'u_seri'
4144	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
4145	varname = 'v_seri'
4146	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
4147	varname = 'q_seri'
4148	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
4149
4150	!================================================================
4151	! Conversion of kinetic and potential energy into heat, for
4152	! parameterisation of subgrid-scale motions
4153	!================================================================
4154
4155	d_t_ec(:,:)=0.
4156	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4157	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
4158	& u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4159	& zmasse,exner,d_t_ec)
4160	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4161
4162	fname = 'after ener_conserv'
4163	varname = 'paprs'
4164	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
4165	varname = 'pplay'
4166	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
4167	varname = 'pphi'
4168	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
4169	varname = 't_seri'
4170	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
4171	varname = 'u_seri'
4172	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
4173	varname = 'v_seri'
4174	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
4175	varname = 'q_seri'
4176	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
4177	varname = 'd_t_ec'
4178	CALL check_var3D(fname, varname, d_t_ec, klon, klev, largest, .FALSE.)
4179
4180	!IM
4181	IF (ip_ebil_phy.ge.1) THEN
4182	ztit='after physic'
4183	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
4184	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
4185	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
4186	!C Comme les tendances de la physique sont ajoute dans la dynamique,
4187	!C on devrait avoir que la variation d'entalpie par la dynamique
4188	!C est egale a la variation de la physique au pas de temps precedent.
4189	!C Donc la somme de ces 2 variations devrait etre nulle.
4190
4191	PRINT *,' LLuis sending: d_h_vcol: ',d_h_vcol,' d_qt ',d_qt,' d_ec ',d_ec
4192
4193	call diagphy(airephy,ztit,ip_ebil_phy &
4194	& , topsw, toplw, solsw, sollw, sens &
4195	& , evap, rain_fall, snow_fall, ztsol &
4196	& , d_h_vcol, d_qt, d_ec &
4197	& , fs_bound, fq_bound )
4198	!C
4199	d_h_vcol_phy=d_h_vcol
4200	!C
4201	END IF
4202	PRINT *,' Lluis 11 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
4203	fname = 'After everything, writting'
4204	varname = 'paprs'
4205	CALL check_var3D(fname, varname, paprs, klon, klev+1, largest, .FALSE.)
4206	varname = 'pplay'
4207	CALL check_var3D(fname, varname, pplay, klon, klev, largest, .FALSE.)
4208	varname = 'pphi'
4209	CALL check_var3D(fname, varname, pphi, klon, klev, largest, .FALSE.)
4210	varname = 't_seri'
4211	CALL check_var3D(fname, varname, t_seri, klon, klev, largest, .FALSE.)
4212	varname = 'u_seri'
4213	CALL check_var3D(fname, varname, u_seri, klon, klev, largest, .FALSE.)
4214	varname = 'v_seri'
4215	CALL check_var3D(fname, varname, v_seri, klon, klev, largest, .FALSE.)
4216	varname = 'q_seri'
4217	CALL check_var3D(fname, varname, q_seri, klon, klev, largest, .FALSE.)
4218
4219	PRINT *,'Lluis Reaching the SORTIES point'
4220
4221	!C
4222	!c=======================================================================
4223	!c SORTIES
4224	!c=======================================================================
4225
4226	!IM Interpolation sur les niveaux de pression du NMC
4227	!c -------------------------------------------------
4228	!c
4229	#include "calcul_STDlev.h"
4230	!c
4231	!c slp sea level pressure
4232	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
4233	!c
4234	!ccc prw = eau precipitable
4235	DO i = 1, klon
4236	prw(i) = 0.
4237	DO k = 1, klev
4238	prw(i) = prw(i) + &
4239	& q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
4240	ENDDO
4241	ENDDO
4242	!c
4243	!IM initialisation + calculs divers diag AMIP2
4244	!c
4245	#include "calcul_divers.h"
4246	!c
4247	IF (type_trac == 'inca') THEN
4248	#ifdef INCA
4249	CALL VTe(VTphysiq)
4250	CALL VTb(VTinca)
4251
4252	CALL chemhook_end ( &
4253	& dtime, &
4254	& pplay, &
4255	& t_seri, &
4256	& tr_seri, &
4257	& nbtr, &
4258	& paprs, &
4259	& q_seri, &
4260	& airephy, &
4261	& pphi, &
4262	& pphis, &
4263	& zx_rh)
4264
4265	CALL VTe(VTinca)
4266	CALL VTb(VTphysiq)
4267	#endif
4268	END IF
4269
4270
4271	!c
4272	!c Convertir les incrementations en tendances
4273	!c
4274	IF (prt_level .GE.10) THEN
4275	print *,'Convertir les incrementations en tendances '
4276	ENDIF
4277	!c
4278	if (mydebug) then
4279	call writefield_phy('u_seri',u_seri,llm)
4280	call writefield_phy('v_seri',v_seri,llm)
4281	call writefield_phy('t_seri',t_seri,llm)
4282	call writefield_phy('q_seri',q_seri,llm)
4283	endif
4284
4285	DO k = 1, klev
4286	DO i = 1, klon
4287	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4288	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4289	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4290	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4291	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4292	ENDDO
4293	ENDDO
4294	!c
4295	IF (nqtot.GE.3) THEN
4296	DO iq = 3, nqtot
4297	DO k = 1, klev
4298	DO i = 1, klon
4299	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4300	ENDDO
4301	ENDDO
4302	ENDDO
4303	ENDIF
4304	!c
4305	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4306	!IM global posePB#include "write_bilKP_ins.h"
4307	!IM global posePB#include "write_bilKP_ave.h"
4308	!c
4309
4310	!c Sauvegarder les valeurs de t et q a la fin de la physique:
4311	!c
4312	DO k = 1, klev
4313	DO i = 1, klon
4314	u_ancien(i,k) = u_seri(i,k)
4315	v_ancien(i,k) = v_seri(i,k)
4316	t_ancien(i,k) = t_seri(i,k)
4317	q_ancien(i,k) = q_seri(i,k)
4318	ENDDO
4319	ENDDO
4320
4321	!!! RomP >>>
4322	IF (nqtot.GE.3) THEN
4323	DO iq = 3, nqtot
4324	DO k = 1, klev
4325	DO i = 1, klon
4326	tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
4327	ENDDO
4328	ENDDO
4329	ENDDO
4330	ENDIF
4331	!!! RomP <<<
4332	!==========================================================================
4333	! Sorties des tendances pour un point particulier
4334	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4335	! pour le debug
4336	! La valeur de igout est attribuee plus haut dans le programme
4337	!==========================================================================
4338
4339	if (prt_level.ge.1) then
4340	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4341	write(lunout,*) &
4342	& 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4343	write(lunout,*) &
4344	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4345	& pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4346	& pctsrf(igout,is_sic)
4347	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4348	do k=1,klev
4349	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4350	& d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4351	& d_t_eva(igout,k)
4352	enddo
4353	write(lunout,*) 'cool,heat'
4354	do k=1,klev
4355	write(lunout,*) cool(igout,k),heat(igout,k)
4356	enddo
4357
4358	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4359	do k=1,klev
4360	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4361	& d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4362	enddo
4363
4364	write(lunout,*) 'd_ps ',d_ps(igout)
4365	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4366	do k=1,klev,klon
4367	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4368	& d_qx(igout,k,1),d_qx(igout,k,2)
4369	enddo
4370	endif
4371
4372	!==========================================================================
4373
4374	!c============================================================
4375	!c Calcul de la temperature potentielle
4376	!c============================================================
4377	DO k = 1, klev
4378	DO i = 1, klon
4379	!cJYG/IM theta en debut du pas de temps
4380	!cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4381	!cJYG/IM theta en fin de pas de temps de physique
4382	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4383	!c thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
4384	!c fth_fonctions.F90 et parkind1.F90
4385	!c sinon thetal=theta
4386	!c thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4387	!c : ql_seri(i,k))
4388	thetal(i,k)=theta(i,k)
4389	ENDDO
4390	ENDDO
4391	!c
4392
4393	!c 22.03.04 BEG
4394	!c=============================================================
4395	!c Ecriture des sorties
4396	!c=============================================================
4397	#ifdef CPP_IOIPSL
4398
4399	!c Recupere des varibles calcule dans differents modules
4400	!c pour ecriture dans histxxx.nc
4401
4402	! Get some variables from module fonte_neige_mod
4403	!L. Fita, LMD. November 2013
4404	!! It is not working after removing starphy.nc and limit.nc?
4405	CALL fonte_neige_get_vars(pctsrf, &
4406	& zxfqcalving, zxfqfonte, zxffonte)
4407
4408
4409
4410	!c=============================================================
4411	! Separation entre thermiques et non thermiques dans les sorties
4412	! de fisrtilp
4413	!c=============================================================
4414
4415	if (iflag_thermals>=1) then
4416	d_t_lscth=0.
4417	d_t_lscst=0.
4418	d_q_lscth=0.
4419	d_q_lscst=0.
4420	do k=1,klev
4421	do i=1,klon
4422	if (ptconvth(i,k)) then
4423	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4424	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4425	else
4426	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4427	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4428	endif
4429	enddo
4430	enddo
4431
4432	do i=1,klon
4433	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4434	plul_th(i)=prfl(i,1)+psfl(i,1)
4435	enddo
4436	endif
4437
4438	PRINT *,' Lluis WRITING outputs! itap: ',itap,' itau_phy: ',itau_phy, &
4439	' itau_w: ',itau_w
4440	PRINT *,' Lluis writting outputs qsol: ',qsol(llp), &
4441	' ftsol: ',ftsol(llp,:)
4442
4443	#include "phys_output_write_new.h"
4444
4445
4446	#ifdef histISCCP
4447	#include "write_histISCCP.h"
4448	#endif
4449
4450	PRINT *,' Lluis WRITING histfiles!'
4451
4452	#ifdef histNMC
4453	#include "write_histhfNMC.h"
4454	#include "write_histdayNMC.h"
4455	#include "write_histmthNMC.h"
4456	#endif
4457
4458	#include "write_histday_seri.h"
4459
4460	#include "write_paramLMDZ_phy.h"
4461
4462	#endif
4463
4464	!c 22.03.04 END
4465	!c
4466	!c====================================================================
4467	!c Si c'est la fin, il faut conserver l'etat de redemarrage
4468	!c====================================================================
4469	!c
4470
4471	!c -----------------------------------------------------------------
4472	!c WSTATS: Saving statistics
4473	!c -----------------------------------------------------------------
4474	!c ("stats" stores and accumulates 8 key variables in file "stats.nc"
4475	!c which can later be used to make the statistic files of the run:
4476	!c "stats") only possible in 3D runs !
4477
4478
4479	IF (callstats) THEN
4480
4481	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
4482	& ,2,paprs(:,1))
4483	call wstats(klon,o_tsol%name,"Surface temperature","K", &
4484	& 2,zxtsol)
4485	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4486	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
4487	& "kg/(s*m2)",2,zx_tmp_fi2d)
4488	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4489	call wstats(klon,o_plul%name,"Large-scale Precip", &
4490	& "kg/(s*m2)",2,zx_tmp_fi2d)
4491	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4492	call wstats(klon,o_pluc%name,"Convective Precip", &
4493	& "kg/(s*m2)",2,zx_tmp_fi2d)
4494	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
4495	& "W/m2",2,solsw)
4496	call wstats(klon,o_soll%name,"IR rad. at surf.", &
4497	& "W/m2",2,sollw)
4498	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4499	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
4500	& "W/m2",2,zx_tmp_fi2d)
4501
4502
4503
4504	call wstats(klon,o_temp%name,"Air temperature","K", &
4505	& 3,t_seri)
4506	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
4507	& 3,u_seri)
4508	call wstats(klon,o_vitv%name,"Meridional wind", &
4509	& "m.s-1",3,v_seri)
4510	call wstats(klon,o_vitw%name,"Vertical wind", &
4511	& "m.s-1",3,omega)
4512	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
4513	& 3,q_seri)
4514
4515
4516
4517	IF(lafin) THEN
4518	write (,) "Writing stats..."
4519	call mkstats(ierr)
4520	ENDIF
4521
4522	ENDIF !if callstats
4523
4524	IF (lafin) THEN
4525	itau_phy = itau_phy + itap
4526	CALL phyredem ("restartphy.nc")
4527	! open(97,form="unformatted",file="finbin")
4528	! write(97) u_seri,v_seri,t_seri,q_seri
4529	! close(97)
4530	!$OMP MASTER
4531	if (read_climoz >= 1) then
4532	if (is_mpi_root) then
4533	call nf95_close(ncid_climoz)
4534	end if
4535	deallocate(press_climoz) ! pointer
4536	end if
4537	!$OMP END MASTER
4538	ENDIF
4539
4540	! first=.false.
4541
4542	! Lluis
4543	PRINT *,' Lluis: ',klev,' UBOUNDS: ',UBOUND(t_seri), UBOUND(u_seri), &
4544	UBOUND(d_q_con), UBOUND(d_t_con)
4545	PRINT *,' Lluis llp ',llp,' itap: ',itap,' zlev t_seri u_seri d_q_con d_t_con_____'
4546	DO i=1,klev
4547	PRINT *,i,t_seri(llp,i), u_seri(llp,i), d_q_con(llp,i), d_t_con(llp,i)
4548	END DO
4549
4550
4551	RETURN
4552	END SUBROUTINE physiq
4553
4554	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4555	IMPLICIT none
4556	!c
4557	!c Calculer et imprimer l'eau totale. A utiliser pour verifier
4558	!c la conservation de l'eau
4559	!c
4560	#include "YOMCST.h"
4561	INTEGER klon,klev
4562	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4563	REAL aire(klon)
4564	REAL qtotal, zx, qcheck
4565	INTEGER i, k
4566	!c
4567	zx = 0.0
4568	DO i = 1, klon
4569	zx = zx + aire(i)
4570	ENDDO
4571	qtotal = 0.0
4572	DO k = 1, klev
4573	DO i = 1, klon
4574	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4575	& *(paprs(i,k)-paprs(i,k+1))/RG
4576	ENDDO
4577	ENDDO
4578	!c
4579	qcheck = qtotal/zx
4580	!c
4581	RETURN
4582	END
4583
4584	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4585	IMPLICIT none
4586	!c
4587	!c Tranformer une variable de la grille physique a
4588	!c la grille d'ecriture
4589	!c
4590	INTEGER nfield,nlon,iim,jjmp1, jjm
4591	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4592	!c
4593	INTEGER i, n, ig
4594	!c
4595	jjm = jjmp1 - 1
4596	DO n = 1, nfield
4597	DO i=1,iim
4598	ecrit(i,n) = fi(1,n)
4599	ecrit(i+jjm*iim,n) = fi(nlon,n)
4600	ENDDO
4601	DO ig = 1, nlon - 2
4602	ecrit(iim+ig,n) = fi(1+ig,n)
4603	ENDDO
4604	ENDDO
4605	RETURN
4606	END SUBROUTINE gr_fi_ecrit
4607
4608	SUBROUTINE check_var(funcn, varn, var, sizev, bigvalue, stoprun)
4609	! Subroutine to check the consistency of a variable
4610	! * NaN value: by definition is variable /= variable
4611	! * bigvalue: threshold for the variable
4612
4613	IMPLICIT NONE
4614
4615	#include "dimensions.h"
4616
4617	INTEGER, INTENT(IN) :: sizev
4618	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
4619	REAL, DIMENSION(sizev), INTENT(IN) :: var
4620	REAL, INTENT(IN) :: bigvalue
4621	LOGICAL, INTENT(IN) :: stoprun
4622
4623	! Local
4624	INTEGER :: i, wrongi, xpt, ypt
4625	CHARACTER(LEN=50) :: errmsg
4626	LOGICAL :: found
4627	REAL, DIMENSION(sizev) :: wrongvalues
4628	INTEGER, DIMENSION(sizev) :: wronggridpt
4629
4630	!!!!!!! Variables
4631	! funcn: at which functino of part of the program variable is checked
4632	! varn: name of the variable
4633	! var: variable to check
4634	! sizev: size of the variable
4635	! bigvalue: biggest attenaible value for the variable
4636	! stoprun: Should the run stop if it founds a problem?
4637
4638	errmsg = 'ERROR -- error -- ERROR -- error'
4639
4640	found = .FALSE.
4641	wrongi = 0
4642	DO i=1,sizev
4643	IF (var(i) /= var(i) .OR. ABS(var(i)) > bigvalue ) THEN
4644	IF (wrongi == 0) found = .TRUE.
4645	wrongi = wrongi + 1
4646	wrongvalues(wrongi) = var(i)
4647	wronggridpt(wrongi) = i
4648	END IF
4649	END DO
4650
4651	IF (found) THEN
4652	PRINT *,TRIM(errmsg)
4653	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
4654	"' is wrong in Nvalues= ",wrongi,' at i (x, y) value___'
4655	DO i=1,wrongi
4656	ypt = INT(wronggridpt(i)/wiim) + 1
4657	xpt = wronggridpt(i) - (ypt-1)*wiim
4658	PRINT *,wronggridpt(i), '(',xpt,', ',ypt,')', wrongvalues(i)
4659	END DO
4660	IF (stoprun) THEN
4661	STOP
4662	END IF
4663	END IF
4664
4665	RETURN
4666
4667	END SUBROUTINE check_var
4668
4669	SUBROUTINE check_var3D(funcn, varn, var, sizev, zsize, bigvalue, stoprun)
4670	! Subroutine to check the consistency of a 3D LMDSZ - variable (klon, klev) !
4671	! * NaN value: by definition is variable /= variable
4672	! * bigvalue: threshold for the variable
4673
4674	IMPLICIT NONE
4675
4676	#include "dimensions.h"
4677
4678	INTEGER, INTENT(IN) :: sizev, zsize
4679	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
4680	REAL, DIMENSION(sizev,zsize), INTENT(IN) :: var
4681	REAL, INTENT(IN) :: bigvalue
4682	LOGICAL, INTENT(IN) :: stoprun
4683
4684	! Local
4685	INTEGER :: i, k, wrongi, xpt, ypt
4686	CHARACTER(LEN=50) :: errmsg
4687	LOGICAL :: found
4688	REAL, DIMENSION(sizev*zsize) :: wrongvalues
4689	INTEGER, DIMENSION(sizev*zsize,2) :: wronggridpt
4690
4691	!!!!!!! Variables
4692	! funcn: at which functino of part of the program variable is checked
4693	! varn: name of the variable
4694	! var: variable to check
4695	! sizev: size of the variable
4696	! zsize: vertical size of the variable
4697	! bigvalue: biggest attenaible value for the variable
4698	! stoprun: Should the run stop if it founds a problem?
4699
4700	errmsg = 'ERROR -- error -- ERROR -- error'
4701
4702	found = .FALSE.
4703	wrongi = 0
4704	DO i=1,sizev
4705	DO k=1,zsize
4706	IF (var(i,k) /= var(i,k) .OR. ABS(var(i,k)) > bigvalue ) THEN
4707	IF (wrongi == 0) found = .TRUE.
4708	wrongi = wrongi + 1
4709	wrongvalues(wrongi) = var(i,k)
4710	wronggridpt(wrongi,1) = i
4711	wronggridpt(wrongi,2) = k
4712	END IF
4713	END DO
4714	END DO
4715
4716	IF (found) THEN
4717	PRINT *,TRIM(errmsg)
4718	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
4719	"' is wrong in Nvalues= ",wrongi,' at i (x,y) k value___'
4720	DO i=1,wrongi
4721	ypt = INT(wronggridpt(i,1)/wiim) + 1
4722	xpt = wronggridpt(i,1) - (ypt-1)*wiim
4723	PRINT *,wronggridpt(i,1), '(',xpt,', ',ypt,')', wronggridpt(i,2), wrongvalues(i)
4724	END DO
4725	IF (stoprun) THEN
4726	STOP
4727	END IF
4728	END IF
4729
4730	RETURN
4731
4732	END SUBROUTINE check_var3D
4733

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

Download in other formats:

Original Format