Context Navigation

physiq.F90 @ 161

Last change on this file since 161 was 161, checked in by lfita, 10 years ago
Fixing REAL/INTEGER issue on calling threscheck_var3D thres is real!
File size: 185.9 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
1304	! L. Fita, LMD August 2014
1305	CHARACTER(LEN=50) :: errmsg, fname, varnamechk
1306	CHARACTER(LEN=50) :: jDS, jHS
1307	LOGICAL :: found
1308	REAL :: largest
1309
1310	llp = 644
1311
1312	fname = 'diagphy'
1313	errmsg = 'ERROR -- error -- ERROR -- error'
1314	largest = 10.e7
1315
1316	!c======================================================================
1317	! Gestion calendrier : mise a jour du module phys_cal_mod
1318	!
1319	CALL phys_cal_update(jD_cur,jH_cur)
1320	WRITE(jDS,'(f10.5)')jD_cur
1321	WRITE(jHS,'(f10.5)')jH_cur
1322
1323	fname = 'Entering in physic at day= '//TRIM(jDS)//' hour: '//TRIM(jHS)
1324	varnamechk = 'paprs'
1325	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
1326	varnamechk = 'pplay'
1327	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
1328	varnamechk = 'pphi'
1329	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
1330	varnamechk = 't'
1331	CALL check_var3D(fname, varnamechk, t, klon, klev, largest, .FALSE.)
1332	varnamechk = 'u'
1333	CALL check_var3D(fname, varnamechk, u, klon, klev, largest, .FALSE.)
1334	varnamechk = 'v'
1335	CALL check_var3D(fname, varnamechk, v, klon, klev, largest, .FALSE.)
1336	varnamechk = 'qv'
1337	CALL check_var3D(fname, varnamechk, qx(:,:,1), klon, klev, largest, .FALSE.)
1338	varnamechk = 'ql'
1339	CALL check_var3D(fname, varnamechk, qx(:,:,2), klon, klev, largest, .FALSE.)
1340	varnamechk = 'qv'
1341	CALL threscheck_var3D(fname, varnamechk, qx(:,:,1), klon, klev, 0.,-1, .FALSE.)
1342
1343
1344	!c======================================================================
1345	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1346	! Utilise notamment en 1D mais peut etre active egalement en 3D
1347	! en imposant la valeur de igout.
1348	!c======================================================================
1349
1350	if (prt_level.ge.1) then
1351	! Lluis
1352	! igout=klon/2+1/klon
1353	igout=llp
1354	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1355	write(lunout,*) &
1356	& 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1357	write(lunout,*) &
1358	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1359
1360	write(lunout,*) 'paprs, play, phi, u, v, t'
1361	do k=1,klev
1362	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1363	& u(igout,k),v(igout,k),t(igout,k)
1364	enddo
1365	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1366	do k=1,klev
1367	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1368	enddo
1369	endif
1370
1371	!c======================================================================
1372
1373	!cym => necessaire pour iflag_con != 2
1374	pmfd(:,:) = 0.
1375	pen_u(:,:) = 0.
1376	pen_d(:,:) = 0.
1377	pde_d(:,:) = 0.
1378	pde_u(:,:) = 0.
1379	aam=0.
1380
1381	torsfc=0.
1382
1383	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1384
1385	if (first) then
1386
1387	!cCR:nvelles variables convection/poches froides
1388
1389	print*, '================================================='
1390	print*, 'Allocation des variables locales et sauvegardees'
1391
1392	call phys_local_var_init
1393
1394	!c
1395	pasphys=pdtphys
1396	!c appel a la lecture du run.def physique
1397	call conf_phys(ok_journe, ok_mensuel, &
1398	& ok_instan, ok_hf, &
1399	& ok_LES, &
1400	& callstats, &
1401	& solarlong0,seuil_inversion, &
1402	& fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1403	& iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1404	& ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
1405	& flag_aerosol, flag_aerosol_strat, new_aod, &
1406	& bl95_b0, bl95_b1, &
1407	!c nv flags pour la convection et les poches froides &
1408	& read_climoz, &
1409	& alp_offset)
1410
1411
1412	!! call phys_state_var_init(read_climoz)
1413	call phys_output_var_init
1414
1415
1416	!L. Fita, LMD. November 2013
1417	! Initializing
1418	! Variables have been already initialilzed from the WRF variables.
1419	! Using on purpose subroutines from physiq_limit_variables_mod
1420	!! CALL neige_initialize()
1421	!! CALL limit_initialize()
1422	!! CALL vars_limit_init(klon)
1423
1424	! L. Fita, LMD January 2014. Initializing output variables
1425	CALL init_ovars_lmdz_NOmodule(klon,klev,nbsrf)
1426
1427	! Allocating restart variables
1428	! ALLOCATE(qsol(klon))
1429	! ALLOCATE(fder(klon))
1430	! ALLOCATE(snow(klon, nbsrf))
1431	! ALLOCATE(qsurf(klon, nbsrf))
1432	! ALLOCATE(evap(klon, nbsrf))
1433	! ALLOCATE(rugos(klon, nbsrf))
1434	! ALLOCATE(agesno(klon, nbsrf))
1435	! ALLOCATE(ftsoil(klon, nsoilmx, nbsrf))
1436	!! ALLOCATE(restart_runoff(klon))
1437
1438	! CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst, &
1439	! rugos_rst, agesno_rst, ftsoil_rst)
1440	! CALL fonte_neige_init(restart_runoff)
1441	! qsol=qsol_rst
1442
1443	print*, '================================================='
1444	!c
1445	dnwd0=0.0
1446	ftd=0.0
1447	fqd=0.0
1448	cin=0.
1449	!cym Attention pbase pas initialise dans concvl !!!!
1450	pbase=0
1451	!IM 180608
1452
1453	itau_con=0
1454	first=.false.
1455
1456	endif ! first
1457
1458	modname = 'physiq'
1459	!IM
1460	IF (ip_ebil_phy.ge.1) THEN
1461	DO i=1,klon
1462	zero_v(i)=0.
1463	END DO
1464	END IF
1465	ok_sync=.TRUE.
1466
1467	IF (debut) THEN
1468	CALL suphel ! initialiser constantes et parametres phys.
1469	ENDIF
1470
1471	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1472
1473
1474	!c======================================================================
1475	! Gestion calendrier : mise a jour du module phys_cal_mod
1476	!
1477	!c CALL phys_cal_update(jD_cur,jH_cur)
1478
1479	!c
1480	!c Si c'est le debut, il faut initialiser plusieurs choses
1481	!c ********
1482	!c
1483
1484	IF (debut) THEN
1485	!rv
1486	!cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1487	!cde la convection a partir des caracteristiques du thermique
1488	wght_th(:,:)=1.
1489	lalim_conv(:)=1
1490	!cRC
1491	ustar(:,:)=0.
1492	u10m(:,:)=0.
1493	v10m(:,:)=0.
1494	rain_con(:)=0.
1495	snow_con(:)=0.
1496	topswai(:)=0.
1497	topswad(:)=0.
1498	solswai(:)=0.
1499	solswad(:)=0.
1500
1501	wmax_th(:)=0.
1502	tau_overturning_th(:)=0.
1503
1504	IF (type_trac == 'inca') THEN
1505	! jg : initialisation jusqu'au ces variables sont dans restart
1506	ccm(:,:,:) = 0.
1507	tau_aero(:,:,:,:) = 0.
1508	piz_aero(:,:,:,:) = 0.
1509	cg_aero(:,:,:,:) = 0.
1510	END IF
1511
1512	rnebcon0(:,:) = 0.0
1513	clwcon0(:,:) = 0.0
1514	rnebcon(:,:) = 0.0
1515	clwcon(:,:) = 0.0
1516
1517	!IM
1518	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1519	!c
1520	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1521	& iflag_coupl,iflag_clos,iflag_wake
1522	print*,'CYCLE_DIURNE', cycle_diurne
1523	!c
1524	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1525	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1526	CALL abort_gcm (modname,abort_message,1)
1527	ENDIF
1528	!c
1529	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1530	abort_message = 'ISCCP-like outputs may be available for KE &
1531	&(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1532	CALL abort_gcm (modname,abort_message,1)
1533	ENDIF
1534	!c
1535	!c Initialiser les compteurs:
1536	!c
1537	itap = 0
1538	itaprad = 0
1539
1540	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1541	!! Un petit travail \`a faire ici.
1542	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1543
1544	if (iflag_pbl>1) then
1545	PRINT*, "Using method MELLOR&YAMADA"
1546	endif
1547
1548	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1549	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1550	! dyn3d
1551	! Attention : la version precedente n'etait pas tres propre.
1552	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1553	! pour obtenir le meme resultat.
1554	dtime=pdtphys
1555	radpas = NINT( 86400./dtime/nbapp_rad)
1556	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1557
1558	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1559
1560	IF (klon_glo==1) THEN
1561	coefh=0. ; coefm=0. ; pbl_tke=0.
1562	coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
1563	PRINT*,'FH WARNING : lignes a supprimer'
1564	ENDIF
1565	!IM begin
1566	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1567	&,ratqs(1,1)
1568	!IM end
1569
1570
1571
1572	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1573	!c
1574	!C on remet le calendrier a zero
1575	!c
1576	IF (raz_date .eq. 1) THEN
1577	itau_phy = 0
1578	ENDIF
1579
1580	!IM cf. AM 081204 BEG
1581	PRINT*,'cycle_diurne3 =',cycle_diurne
1582	!IM cf. AM 081204 END
1583	!c
1584	CALL printflag( tabcntr0,radpas,ok_journe, &
1585	& ok_instan, ok_region )
1586	!c
1587	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1588	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1589	& pdtphys
1590	abort_message='Pas physique n est pas correct '
1591	! call abort_gcm(modname,abort_message,1)
1592	dtime=pdtphys
1593	ENDIF
1594	IF (nlon .NE. klon) THEN
1595	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1596	& klon
1597	abort_message='nlon et klon ne sont pas coherents'
1598	call abort_gcm(modname,abort_message,1)
1599	ENDIF
1600	IF (nlev .NE. klev) THEN
1601	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1602	& klev
1603	abort_message='nlev et klev ne sont pas coherents'
1604	call abort_gcm(modname,abort_message,1)
1605	ENDIF
1606	!c
1607	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1608	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1609	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1610	abort_message='Nbre d appels au rayonnement insuffisant'
1611	call abort_gcm(modname,abort_message,1)
1612	ENDIF
1613	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1614	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1615	& ok_cvl
1616	!c
1617	!cKE43
1618	!c Initialisation pour la convection de K.E. (sb):
1619	IF (iflag_con.GE.3) THEN
1620
1621	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1622	WRITE(lunout,*) &
1623	& "On va utiliser le melange convectif des traceurs qui"
1624	WRITE(lunout,*)"est calcule dans convect4.3"
1625	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1626
1627	DO i = 1, klon
1628	ema_cbmf(i) = 0.
1629	ema_pcb(i) = 0.
1630	ema_pct(i) = 0.
1631	!c ema_workcbmf(i) = 0.
1632	ENDDO
1633	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1634	DO i = 1, klon
1635	ibas_con(i) = 1
1636	itop_con(i) = 1
1637	ENDDO
1638	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1639	!c===============================================================================
1640	!cCR:04.12.07: initialisations poches froides
1641	!c Controle de ALE et ALP pour la fermeture convective (jyg)
1642	if (iflag_wake>=1) then
1643	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1644	& ,alp_bl_prescr, ale_bl_prescr)
1645	!c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1646	!c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1647	endif
1648
1649	do i = 1,klon
1650	Ale_bl(i)=0.
1651	Alp_bl(i)=0.
1652	enddo
1653
1654	!c================================================================================
1655	!IM stations CFMIP
1656	nCFMIP=npCFMIP
1657	OPEN(98,file='npCFMIP_param.data',status='old', &
1658	& form='formatted',iostat=iostat)
1659	if (iostat == 0) then
1660	READ(98,*,end=998) nCFMIP
1661	998 CONTINUE
1662	CLOSE(98)
1663	CONTINUE
1664	IF(nCFMIP.GT.npCFMIP) THEN
1665	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1666	CALL abort
1667	else
1668	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1669	ENDIF
1670
1671	!c
1672	ALLOCATE(tabCFMIP(nCFMIP))
1673	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1674	ALLOCATE(tabijGCM(nCFMIP))
1675	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1676	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1677	!c
1678	!c lecture des nCFMIP stations CFMIP, de leur numero
1679	!c et des coordonnees geographiques lonCFMIP, latCFMIP
1680	!c
1681	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1682	&lonCFMIP, latCFMIP)
1683	!c
1684	!c identification des
1685	!c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
1686	!c 2) indices points tabijGCM de la grille physique 1d sur klon points
1687	!c 3) indices iGCM, jGCM de la grille physique 2d
1688	!c
1689	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1690	&tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1691	!c
1692	else
1693	ALLOCATE(tabijGCM(0))
1694	ALLOCATE(lonGCM(0), latGCM(0))
1695	ALLOCATE(iGCM(0), jGCM(0))
1696	end if
1697	else
1698	ALLOCATE(tabijGCM(0))
1699	ALLOCATE(lonGCM(0), latGCM(0))
1700	ALLOCATE(iGCM(0), jGCM(0))
1701	ENDIF
1702
1703	DO i=1,klon
1704	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1705	ENDDO
1706
1707	!c34EK
1708	IF (ok_orodr) THEN
1709
1710	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1711	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1712	! justement quand ok_orodr = false.
1713	! ce rugoro est utilise par la couche limite et fait double emploi
1714	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1715	! DO i=1,klon
1716	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1717	! ENDDO
1718	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1719	IF (ok_strato) THEN
1720	CALL SUGWD_strato(klon,klev,paprs,pplay)
1721	ELSE
1722	CALL SUGWD(klon,klev,paprs,pplay)
1723	ENDIF
1724
1725	DO i=1,klon
1726	zuthe(i)=0.
1727	zvthe(i)=0.
1728	if(zstd(i).gt.10.)then
1729	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1730	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1731	endif
1732	ENDDO
1733	ENDIF
1734	!c
1735	!c
1736	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1737	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1738	& lmt_pas
1739	!c
1740	capemaxcels = 't_max(X)'
1741	t2mincels = 't_min(X)'
1742	t2maxcels = 't_max(X)'
1743	tinst = 'inst(X)'
1744	tave = 'ave(X)'
1745	!IM cf. AM 081204 BEG
1746	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1747	!IM cf. AM 081204 END
1748	!c
1749	!c=============================================================
1750	!c Initialisation des sorties
1751	!c=============================================================
1752
1753	#ifdef CPP_IOIPSL
1754
1755	!$OMP MASTER
1756	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
1757	& iGCM,jGCM,lonGCM,latGCM, &
1758	& jjmp1,nlevSTD,clevSTD, &
1759	& nbteta, ctetaSTD, dtime,ok_veget, &
1760	& type_ocean,iflag_pbl,ok_mensuel,ok_journe, &
1761	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1762	& read_climoz, phys_out_filestations, &
1763	& new_aod, aerosol_couple, &
1764	& flag_aerosol_strat )
1765	!$OMP END MASTER
1766	!$OMP BARRIER
1767
1768	#ifdef histISCCP
1769	#include "ini_histISCCP.h"
1770	#endif
1771
1772	#ifdef histNMC
1773	#include "ini_histhfNMC.h"
1774	#include "ini_histdayNMC.h"
1775	#include "ini_histmthNMC.h"
1776	#endif
1777
1778	#include "ini_histday_seri.h"
1779
1780	#include "ini_paramLMDZ_phy.h"
1781
1782	#endif
1783	ecrit_reg = ecrit_reg * un_jour
1784	ecrit_tra = ecrit_tra * un_jour
1785
1786	!cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1787	date0 = jD_ref
1788	WRITE(,) 'physiq date0 : ',date0
1789	!c
1790	!c
1791	!c
1792	!c Prescrire l'ozone dans l'atmosphere
1793	!c
1794	!c
1795	!cc DO i = 1, klon
1796	!cc DO k = 1, klev
1797	!cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1798	!cc ENDDO
1799	!cc ENDDO
1800	!c
1801	IF (type_trac == 'inca') THEN
1802	#ifdef INCA
1803	CALL VTe(VTphysiq)
1804	CALL VTb(VTinca)
1805	! iii = MOD(NINT(xjour),360)
1806	! calday = REAL(iii) + jH_cur
1807	calday = REAL(days_elapsed) + jH_cur
1808	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1809
1810	CALL chemini( &
1811	& rg, &
1812	& ra, &
1813	& airephy, &
1814	& rlat, &
1815	& rlon, &
1816	& presnivs, &
1817	& calday, &
1818	& klon, &
1819	& nqtot, &
1820	& pdtphys, &
1821	& annee_ref, &
1822	& day_ref, &
1823	& itau_phy)
1824
1825	CALL VTe(VTinca)
1826	CALL VTb(VTphysiq)
1827	#endif
1828	END IF
1829	!c
1830	!c
1831	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1832	! Nouvelle initialisation pour le rayonnement RRTM
1833	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1834
1835	call iniradia(klon,klev,paprs(1,1:klev+1))
1836
1837	!$OMP single
1838	if (read_climoz >= 1) then
1839	call open_climoz(ncid_climoz, press_climoz)
1840	END IF
1841	!$OMP end single
1842	!c
1843	!IM betaCRF
1844	pfree=70000. !Pa
1845	beta_pbl=1.
1846	beta_free=1.
1847	lon1_beta=-180.
1848	lon2_beta=+180.
1849	lat1_beta=90.
1850	lat2_beta=-90.
1851	mskocean_beta=.FALSE.
1852
1853	OPEN(99,file='beta_crf.data',status='old', &
1854	& form='formatted',err=9999)
1855	READ(99,*,end=9998) pfree
1856	READ(99,*,end=9998) beta_pbl
1857	READ(99,*,end=9998) beta_free
1858	READ(99,*,end=9998) lon1_beta
1859	READ(99,*,end=9998) lon2_beta
1860	READ(99,*,end=9998) lat1_beta
1861	READ(99,*,end=9998) lat2_beta
1862	READ(99,*,end=9998) mskocean_beta
1863	9998 Continue
1864	CLOSE(99)
1865	9999 Continue
1866	WRITE(,)'pfree=',pfree
1867	WRITE(,)'beta_pbl=',beta_pbl
1868	WRITE(,)'beta_free=',beta_free
1869	WRITE(,)'lon1_beta=',lon1_beta
1870	WRITE(,)'lon2_beta=',lon2_beta
1871	WRITE(,)'lat1_beta=',lat1_beta
1872	WRITE(,)'lat2_beta=',lat2_beta
1873	WRITE(,)'mskocean_beta=',mskocean_beta
1874	ENDIF
1875	!
1876	! ************** Fin de IF ( debut ) *************
1877	!
1878	!
1879	! Incrementer le compteur de la physique
1880	!
1881	itap = itap + 1
1882	WRITE(jDS,'(I5)')itap
1883	fname = 'After itap= '//TRIM(jDS)
1884	varnamechk = 'paprs'
1885	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
1886	varnamechk = 'pplay'
1887	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
1888	varnamechk = 'pphi'
1889	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
1890	varnamechk = 't'
1891	CALL check_var3D(fname, varnamechk, t, klon, klev, largest, .FALSE.)
1892	varnamechk = 'u'
1893	CALL check_var3D(fname, varnamechk, u, klon, klev, largest, .FALSE.)
1894	varnamechk = 'v'
1895	CALL check_var3D(fname, varnamechk, v, klon, klev, largest, .FALSE.)
1896	varnamechk = 'qv'
1897	CALL check_var3D(fname, varnamechk, qx(:,:,1), klon, klev, largest, .FALSE.)
1898	varnamechk = 'ql'
1899	CALL check_var3D(fname, varnamechk, qx(:,:,2), klon, klev, largest, .FALSE.)
1900
1901	!c
1902	!
1903	! Update fraction of the sub-surfaces (pctsrf) and
1904	! initialize, where a new fraction has appeared, all variables depending
1905	! on the surface fraction.
1906	!
1907	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1908	& pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1909
1910	! Update time and other variables in Reprobus
1911	IF (type_trac == 'repr') THEN
1912	#ifdef REPROBUS
1913	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1914	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1915	CALL Rtime(debut)
1916	#endif
1917	END IF
1918
1919
1920	! Tendances bidons pour les processus qui n'affectent pas certaines
1921	! variables.
1922	du0(:,:)=0.
1923	dv0(:,:)=0.
1924	dq0(:,:)=0.
1925	dql0(:,:)=0.
1926	!c
1927	!c Mettre a zero des variables de sortie (pour securite)
1928	!c
1929	DO i = 1, klon
1930	d_ps(i) = 0.0
1931	ENDDO
1932	DO k = 1, klev
1933	DO i = 1, klon
1934	d_t(i,k) = 0.0
1935	d_u(i,k) = 0.0
1936	d_v(i,k) = 0.0
1937	ENDDO
1938	ENDDO
1939	DO iq = 1, nqtot
1940	DO k = 1, klev
1941	DO i = 1, klon
1942	d_qx(i,k,iq) = 0.0
1943	ENDDO
1944	ENDDO
1945	ENDDO
1946	da(:,:)=0.
1947	mp(:,:)=0.
1948	phi(:,:,:)=0.
1949	! RomP >>>
1950	phi2(:,:,:)=0.
1951	beta_prec_fisrt(:,:)=0.
1952	beta_prec(:,:)=0.
1953	epmlmMm(:,:,:)=0.
1954	eplaMm(:,:)=0.
1955	d1a(:,:)=0.
1956	dam(:,:)=0.
1957	pmflxr=0.
1958	pmflxs=0.
1959	! RomP <<<
1960
1961	!c
1962	!c Ne pas affecter les valeurs entrees de u, v, h, et q
1963	!c
1964	DO k = 1, klev
1965	DO i = 1, klon
1966	t_seri(i,k) = t(i,k)
1967	u_seri(i,k) = u(i,k)
1968	v_seri(i,k) = v(i,k)
1969	q_seri(i,k) = qx(i,k,ivap)
1970	ql_seri(i,k) = qx(i,k,iliq)
1971	qs_seri(i,k) = 0.
1972	ENDDO
1973	ENDDO
1974	tke0(:,:)=pbl_tke(:,:,is_ave)
1975	IF (nqtot.GE.3) THEN
1976	DO iq = 3, nqtot
1977	DO k = 1, klev
1978	DO i = 1, klon
1979	tr_seri(i,k,iq-2) = qx(i,k,iq)
1980	ENDDO
1981	ENDDO
1982	ENDDO
1983	ELSE
1984	DO k = 1, klev
1985	DO i = 1, klon
1986	tr_seri(i,k,1) = 0.0
1987	ENDDO
1988	ENDDO
1989	ENDIF
1990	!C
1991	DO i = 1, klon
1992	ztsol(i) = 0.
1993	ENDDO
1994	DO nsrf = 1, nbsrf
1995	DO i = 1, klon
1996	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1997	ENDDO
1998	ENDDO
1999	!IM
2000	IF (ip_ebil_phy.ge.1) THEN
2001	ztit='after dynamic'
2002	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
2003	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2004	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2005	!C Comme les tendances de la physique sont ajoute dans la dynamique,
2006	!C on devrait avoir que la variation d'entalpie par la dynamique
2007	!C est egale a la variation de la physique au pas de temps precedent.
2008	!C Donc la somme de ces 2 variations devrait etre nulle.
2009	call diagphy(airephy,ztit,ip_ebil_phy &
2010	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2011	& , zero_v, zero_v, zero_v, ztsol &
2012	& , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
2013	& , fs_bound, fq_bound )
2014	END IF
2015	PRINT *,' Lluis 1 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2016	fname = 'After dynamic'
2017	varnamechk = 'paprs'
2018	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
2019	varnamechk = 'pplay'
2020	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
2021	varnamechk = 'pphi'
2022	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
2023	varnamechk = 't_seri'
2024	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
2025	varnamechk = 'u_seri'
2026	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
2027	varnamechk = 'v_seri'
2028	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
2029	varnamechk = 'q_seri'
2030	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
2031
2032	!c Diagnostiquer la tendance dynamique
2033	!c
2034	IF (ancien_ok) THEN
2035	DO k = 1, klev
2036	DO i = 1, klon
2037	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
2038	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
2039	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
2040	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
2041	ENDDO
2042	ENDDO
2043	!!! RomP >>> td dyn traceur
2044	IF (nqtot.GE.3) THEN
2045	DO iq = 3, nqtot
2046	DO k = 1, klev
2047	DO i = 1, klon
2048	d_tr_dyn(i,k,iq-2)= &
2049	& (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
2050	! iiq=niadv(iq)
2051	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
2052	ENDDO
2053	ENDDO
2054	ENDDO
2055	ENDIF
2056	!!! RomP <<<
2057	ELSE
2058	DO k = 1, klev
2059	DO i = 1, klon
2060	d_u_dyn(i,k) = 0.0
2061	d_v_dyn(i,k) = 0.0
2062	d_t_dyn(i,k) = 0.0
2063	d_q_dyn(i,k) = 0.0
2064	ENDDO
2065	ENDDO
2066	!!! RomP >>> td dyn traceur
2067	IF (nqtot.GE.3) THEN
2068	DO iq = 3, nqtot
2069	DO k = 1, klev
2070	DO i = 1, klon
2071	d_tr_dyn(i,k,iq-2)= 0.0
2072	ENDDO
2073	ENDDO
2074	ENDDO
2075	ENDIF
2076	!!! RomP <<<
2077	ancien_ok = .TRUE.
2078	ENDIF
2079	!c
2080	!c Ajouter le geopotentiel du sol:
2081	!c
2082	DO k = 1, klev
2083	DO i = 1, klon
2084	zphi(i,k) = pphi(i,k) + pphis(i)
2085	ENDDO
2086	ENDDO
2087	!c
2088	!c Verifier les temperatures
2089	!c
2090	!IM BEG
2091	IF (check) THEN
2092	amn=MIN(ftsol(1,is_ter),1000.)
2093	amx=MAX(ftsol(1,is_ter),-1000.)
2094	DO i=2, klon
2095	amn=MIN(ftsol(i,is_ter),amn)
2096	amx=MAX(ftsol(i,is_ter),amx)
2097	ENDDO
2098	!c
2099	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2100	ENDIF !(check) THEN
2101
2102	!IM END
2103	!c
2104
2105	CALL hgardfou(t_seri,ftsol,'debutphy')
2106	!c
2107	!IM BEG
2108	IF (check) THEN
2109	amn=MIN(ftsol(1,is_ter),1000.)
2110	amx=MAX(ftsol(1,is_ter),-1000.)
2111	DO i=2, klon
2112	amn=MIN(ftsol(i,is_ter),amn)
2113	amx=MAX(ftsol(i,is_ter),amx)
2114	ENDDO
2115	!c
2116	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2117	ENDIF !(check) THEN
2118
2119	!IM END
2120	!c
2121	!c Mettre en action les conditions aux limites (albedo, sst, etc.).
2122	!c Prescrire l'ozone et calculer l'albedo sur l'ocean.
2123	!c
2124	if (read_climoz >= 1) then
2125	!C Ozone from a file
2126	! Update required ozone index:
2127	ro3i = int((days_elapsed + jh_cur - jh_1jan) &
2128	& / ioget_year_len(year_cur) * 360.) + 1
2129	if (ro3i == 361) ro3i = 360
2130	!C (This should never occur, except perhaps because of roundup
2131	!C error. See documentation.)
2132	if (ro3i /= co3i) then
2133	!C Update ozone field:
2134	if (read_climoz == 1) then
2135	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
2136	& press_in_edg=press_climoz, paprs=paprs, v3=wo)
2137	else
2138	!C read_climoz == 2
2139	call regr_pr_av(ncid_climoz, &
2140	& (/"tro3 ", "tro3_daylight"/), &
2141	& julien=ro3i, press_in_edg=press_climoz, paprs=paprs, &
2142	& v3=wo)
2143	end if
2144	! Convert from mole fraction of ozone to column density of ozone in a
2145	! cell, in kDU:
2146	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) &
2147	& * rmo3 / rmd * zmasse / dobson_u / 1e3
2148	!C (By regridding ozone values for LMDZ only once every 360th of
2149	!C year, we have already neglected the variation of pressure in one
2150	!C 360th of year. So do not recompute "wo" at each time step even if
2151	!C "zmasse" changes a little.)
2152	co3i = ro3i
2153	end if
2154	elseif (MOD(itap-1,lmt_pas) == 0) THEN
2155	!C Once per day, update ozone from Royer:
2156	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
2157	ENDIF
2158
2159	!c
2160	!c Re-evaporer l'eau liquide nuageuse
2161	!c
2162	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
2163	DO i = 1, klon
2164	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2165	!c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2166	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2167	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
2168	zb = MAX(0.0,ql_seri(i,k))
2169	za = - MAX(0.0,ql_seri(i,k)) &
2170	& * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
2171	t_seri(i,k) = t_seri(i,k) + za
2172	q_seri(i,k) = q_seri(i,k) + zb
2173	ql_seri(i,k) = 0.0
2174	d_t_eva(i,k) = za
2175	d_q_eva(i,k) = zb
2176	ENDDO
2177	ENDDO
2178	!IM
2179	IF (ip_ebil_phy.ge.2) THEN
2180	ztit='after reevap'
2181	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
2182	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2183	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2184	call diagphy(airephy,ztit,ip_ebil_phy &
2185	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2186	& , zero_v, zero_v, zero_v, ztsol &
2187	& , d_h_vcol, d_qt, d_ec &
2188	& , fs_bound, fq_bound )
2189	!C
2190	END IF
2191	PRINT *,' Lluis 2 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2192	fname = 'After reevap'
2193	varnamechk = 'paprs'
2194	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
2195	varnamechk = 'pplay'
2196	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
2197	varnamechk = 'pphi'
2198	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
2199	varnamechk = 't_seri'
2200	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
2201	varnamechk = 'u_seri'
2202	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
2203	varnamechk = 'v_seri'
2204	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
2205	varnamechk = 'q_seri'
2206	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
2207	!c
2208	!c=========================================================================
2209	! Calculs de l'orbite.
2210	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2211	! doit donc etre plac\'e avant radlwsw et pbl_surface
2212
2213	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2214	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2215	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2216	!
2217	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2218	! solarlong0
2219	if (solarlong0<-999.) then
2220	if (new_orbit) then
2221	! calcul selon la routine utilisee pour les planetes
2222	call solarlong(day_since_equinox, zlongi, dist)
2223	else
2224	! calcul selon la routine utilisee pour l'AR4
2225	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2226	endif
2227	else
2228	zlongi=solarlong0 ! longitude solaire vraie
2229	dist=1. ! distance au soleil / moyenne
2230	endif
2231	if(prt_level.ge.1) &
2232	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist,&
2233	' jD_cur: ',jD_cur,' jH_cur: ',jH_cur,' days_elapsed: ', &
2234	REAL(days_elapsed+1)
2235
2236
2237	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2238	! Calcul de l'ensoleillement :
2239	! ============================
2240	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2241	! l'annee a partir d'une formule analytique.
2242	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2243	! non nul aux poles.
2244	IF (abs(solarlong0-1000.)<1.e-4) then
2245	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
2246	ELSE
2247	! Avec ou sans cycle diurne
2248	IF (cycle_diurne) THEN
2249	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
2250	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
2251	ELSE
2252	CALL angle(zlongi, rlat, fract, rmu0)
2253	ENDIF
2254	ENDIF
2255
2256	if (mydebug) then
2257	call writefield_phy('u_seri',u_seri,llm)
2258	call writefield_phy('v_seri',v_seri,llm)
2259	call writefield_phy('t_seri',t_seri,llm)
2260	call writefield_phy('q_seri',q_seri,llm)
2261	endif
2262
2263	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2264	!c Appel au pbl_surface : Planetary Boudary Layer et Surface
2265	!c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2266	!c turbulent du couche limit.
2267	!c
2268	!c Certains varibales de sorties de pbl_surface sont utiliser que pour
2269	!c ecriture des fihiers hist_XXXX.nc, ces sont :
2270	!c qsol, zq2m, s_pblh, s_lcl,
2271	!c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2272	!c s_therm, s_trmb1, s_trmb2, s_trmb3,
2273	!c zxrugs, zu10m, zv10m, fder,
2274	!c zxqsurf, rh2m, zxfluxu, zxfluxv,
2275	!c frugs, agesno, fsollw, fsolsw,
2276	!c d_ts, fevap, fluxlat, t2m,
2277	!c wfbils, wfbilo, fluxt, fluxu, fluxv,
2278	!c
2279	!c Certains ne sont pas utiliser du tout :
2280	!c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2281	!c
2282
2283	!c Calcul de l'humidite de saturation au niveau du sol
2284
2285
2286
2287	if (iflag_pbl/=0) then
2288
2289	CALL pbl_surface( &
2290	& dtime, date0, itap, days_elapsed+1, &
2291	& debut, lafin, &
2292	& rlon, rlat, rugoro, rmu0, &
2293	& rain_fall, snow_fall, solsw, sollw, &
2294	& t_seri, q_seri, u_seri, v_seri, &
2295	& pplay, paprs, pctsrf, &
2296	& ftsol, falb1, falb2, ustar, u10m, v10m, &
2297	& sollwdown, cdragh, cdragm, u1, v1, &
2298	& albsol1, albsol2, sens, evap, &
2299	& zxtsol, zxfluxlat, zt2m, qsat2m, &
2300	& d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2301	& coefh, coefm, slab_wfbils, &
2302	& qsol, zq2m, s_pblh, s_lcl, &
2303	& s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2304	& s_therm, s_trmb1, s_trmb2, s_trmb3, &
2305	& zxrugs, zustar, zu10m, zv10m, fder, &
2306	& zxqsurf, rh2m, zxfluxu, zxfluxv, &
2307	& frugs, agesno, fsollw, fsolsw, &
2308	& d_ts, fevap, fluxlat, t2m, &
2309	& wfbils, wfbilo, fluxt, fluxu, fluxv, &
2310	& dsens, devap, zxsnow, &
2311	& zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
2312
2313	fname = 'After pbl_surface'
2314	varnamechk = 'paprs'
2315	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
2316	varnamechk = 'pplay'
2317	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
2318	varnamechk = 'pphi'
2319	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
2320	varnamechk = 't_seri'
2321	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
2322	varnamechk = 'u_seri'
2323	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
2324	varnamechk = 'v_seri'
2325	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
2326	varnamechk = 'q_seri'
2327	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
2328
2329	!-----------------------------------------------------------------------------------------
2330	! ajout des tendances de la diffusion turbulente
2331	CALL add_phys_tend &
2332	& (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
2333	!-----------------------------------------------------------------------------------------
2334	if (mydebug) then
2335	call writefield_phy('u_seri',u_seri,llm)
2336	call writefield_phy('v_seri',v_seri,llm)
2337	call writefield_phy('t_seri',t_seri,llm)
2338	call writefield_phy('q_seri',q_seri,llm)
2339	endif
2340
2341	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2342	& t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2343
2344
2345	IF (ip_ebil_phy.ge.2) THEN
2346	ztit='after surface_main'
2347	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2348	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2349	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2350	call diagphy(airephy,ztit,ip_ebil_phy &
2351	& , zero_v, zero_v, zero_v, zero_v, sens &
2352	& , evap , zero_v, zero_v, ztsol &
2353	& , d_h_vcol, d_qt, d_ec &
2354	& , fs_bound, fq_bound )
2355	END IF
2356
2357	ENDIF
2358	PRINT *,' Lluis 3 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2359	fname = 'After surface_main'
2360	varnamechk = 'paprs'
2361	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
2362	varnamechk = 'pplay'
2363	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
2364	varnamechk = 'pphi'
2365	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
2366	varnamechk = 't_seri'
2367	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
2368	varnamechk = 'u_seri'
2369	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
2370	varnamechk = 'v_seri'
2371	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
2372	varnamechk = 'q_seri'
2373	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
2374
2375	!c =================================================================== c
2376	!c Calcul de Qsat
2377
2378	DO k = 1, klev
2379	DO i = 1, klon
2380	zx_t = t_seri(i,k)
2381	IF (thermcep) THEN
2382	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2383	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2384	zx_qs = MIN(0.5,zx_qs)
2385	zcor = 1./(1.-retv*zx_qs)
2386	zx_qs = zx_qs*zcor
2387	ELSE
2388	IF (zx_t.LT.t_coup) THEN
2389	zx_qs = qsats(zx_t)/pplay(i,k)
2390	ELSE
2391	zx_qs = qsatl(zx_t)/pplay(i,k)
2392	ENDIF
2393	ENDIF
2394	zqsat(i,k)=zx_qs
2395	ENDDO
2396	ENDDO
2397
2398	if (prt_level.ge.1) then
2399	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2400	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2401	endif
2402	!c
2403	!c Appeler la convection (au choix)
2404	!c
2405	DO k = 1, klev
2406	DO i = 1, klon
2407	conv_q(i,k) = d_q_dyn(i,k) &
2408	& + d_q_vdf(i,k)/dtime
2409	conv_t(i,k) = d_t_dyn(i,k) &
2410	& + d_t_vdf(i,k)/dtime
2411	ENDDO
2412	ENDDO
2413	IF (check) THEN
2414	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2415	WRITE(lunout,*) "avantcon=", za
2416	ENDIF
2417	zx_ajustq = .FALSE.
2418	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2419	IF (zx_ajustq) THEN
2420	DO i = 1, klon
2421	z_avant(i) = 0.0
2422	ENDDO
2423	DO k = 1, klev
2424	DO i = 1, klon
2425	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2426	& *(paprs(i,k)-paprs(i,k+1))/RG
2427	ENDDO
2428	ENDDO
2429	ENDIF
2430
2431	!c Calcule de vitesse verticale a partir de flux de masse verticale
2432	DO k = 1, klev
2433	DO i = 1, klon
2434	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2435	END DO
2436	END DO
2437	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2438	& omega(igout, :)
2439
2440	IF (iflag_con.EQ.1) THEN
2441	abort_message ='reactiver le call conlmd dans physiq.F'
2442	CALL abort_gcm (modname,abort_message,1)
2443	!c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2444	!c . d_t_con, d_q_con,
2445	!c . rain_con, snow_con, ibas_con, itop_con)
2446	ELSE IF (iflag_con.EQ.2) THEN
2447	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
2448	& conv_t, conv_q, -evap, omega, &
2449	& d_t_con, d_q_con, rain_con, snow_con, &
2450	& pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2451	& kcbot, kctop, kdtop, pmflxr, pmflxs)
2452	d_u_con = 0.
2453	d_v_con = 0.
2454
2455	WHERE (rain_con < 0.) rain_con = 0.
2456	WHERE (snow_con < 0.) snow_con = 0.
2457	DO i = 1, klon
2458	ibas_con(i) = klev+1 - kcbot(i)
2459	itop_con(i) = klev+1 - kctop(i)
2460	ENDDO
2461	ELSE IF (iflag_con.GE.3) THEN
2462	!c nb of tracers for the KE convection:
2463	!c MAF la partie traceurs est faite dans phytrac
2464	!c on met ntra=1 pour limiter les appels mais on peut
2465	!c supprimer les calculs / ftra.
2466	ntra = 1
2467
2468	!c=====================================================================================
2469	!cajout pour la parametrisation des poches froides:
2470	!ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2471	do k=1,klev
2472	do i=1,klon
2473	if (iflag_wake>=1) then
2474	t_wake(i,k) = t_seri(i,k) &
2475	& +(1-wake_s(i))*wake_deltat(i,k)
2476	q_wake(i,k) = q_seri(i,k) &
2477	& +(1-wake_s(i))*wake_deltaq(i,k)
2478	t_undi(i,k) = t_seri(i,k) &
2479	& -wake_s(i)*wake_deltat(i,k)
2480	q_undi(i,k) = q_seri(i,k) &
2481	& -wake_s(i)*wake_deltaq(i,k)
2482	else
2483	t_wake(i,k) = t_seri(i,k)
2484	q_wake(i,k) = q_seri(i,k)
2485	t_undi(i,k) = t_seri(i,k)
2486	q_undi(i,k) = q_seri(i,k)
2487	endif
2488	enddo
2489	enddo
2490
2491	!cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2492	!cc-- pour le soulevement des particules dans le modele convectif
2493	!c
2494	do i = 1,klon
2495	ALE(i) = 0.
2496	ALP(i) = 0.
2497	enddo
2498	!c
2499	!ccalcul de ale_wake et alp_wake
2500	if (iflag_wake>=1) then
2501	if (itap .le. it_wape_prescr) then
2502	do i = 1,klon
2503	ale_wake(i) = wape_prescr
2504	alp_wake(i) = fip_prescr
2505	enddo
2506	else
2507	do i = 1,klon
2508	!cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2509	!ccc ale_wake(i) = 0.5wake_cstar(i)*2
2510	ale_wake(i) = wake_pe(i)
2511	alp_wake(i) = wake_fip(i)
2512	enddo
2513	endif
2514	else
2515	do i = 1,klon
2516	ale_wake(i) = 0.
2517	alp_wake(i) = 0.
2518	enddo
2519	endif
2520	!ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2521	!cdans le thermique sinon
2522	if (iflag_coupl.eq.0) then
2523	if (debut.and.prt_level.gt.9) &
2524	& WRITE(lunout,*)'ALE et ALP imposes'
2525	do i = 1,klon
2526	!con ne couple que ale
2527	!c ALE(i) = max(ale_wake(i),Ale_bl(i))
2528	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2529	!con ne couple que alp
2530	!c ALP(i) = alp_wake(i) + Alp_bl(i)
2531	ALP(i) = alp_wake(i) + alp_bl_prescr
2532	enddo
2533	else
2534	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2535	! do i = 1,klon
2536	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2537	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2538	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2539	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2540	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2541	! enddo
2542
2543	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2544	! Modif FH 2010/04/27. Sans doute temporaire.
2545	! Deux options pour le alp_offset : constant si >?? 0 ou proportionnel ??a
2546	! w si <0
2547	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2548	do i = 1,klon
2549	ALE(i) = max(ale_wake(i),Ale_bl(i))
2550	!ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
2551	if (iflag_trig_bl.ge.1) then
2552	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2553	endif
2554	!ccc fin nrlmd le 10/04/2012
2555	if (alp_offset>=0.) then
2556	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2557	else
2558	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2559	if (alp(i)<0.) then
2560	print*,'ALP ',alp(i),alp_wake(i) &
2561	& ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2562	endif
2563	endif
2564	enddo
2565	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2566
2567	endif
2568	do i=1,klon
2569	if (alp(i)>alp_max) then
2570	IF(prt_level>9)WRITE(lunout,*) &
2571	& 'WARNING SUPER ALP (seuil=',alp_max, &
2572	& '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2573	alp(i)=alp_max
2574	endif
2575	if (ale(i)>ale_max) then
2576	IF(prt_level>9)WRITE(lunout,*) &
2577	& 'WARNING SUPER ALE (seuil=',ale_max, &
2578	& '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2579	ale(i)=ale_max
2580	endif
2581	enddo
2582
2583	!cfin calcul ale et alp
2584	!c=================================================================================================
2585
2586
2587	!c sb, oct02:
2588	!c Schema de convection modularise et vectorise:
2589	!c (driver commun aux versions 3 et 4)
2590	!c
2591	IF (ok_cvl) THEN ! new driver for convectL
2592
2593	IF (type_trac == 'repr') THEN
2594	nbtr_tmp=ntra
2595	ELSE
2596	nbtr_tmp=nbtr
2597	END IF
2598	CALL concvl (iflag_con,iflag_clos, &
2599	& dtime,paprs,pplay,t_undi,q_undi, &
2600	& t_wake,q_wake,wake_s, &
2601	& u_seri,v_seri,tr_seri,nbtr_tmp, &
2602	& ALE,ALP, &
2603	& ema_work1,ema_work2, &
2604	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2605	& rain_con, snow_con, ibas_con, itop_con, sigd, &
2606	& ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2607	& Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2608	& pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2609	! RomP >>>
2610	!! . pmflxr,pmflxs,da,phi,mp,
2611	!! . ftd,fqd,lalim_conv,wght_th) &
2612	& pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2613	& ftd,fqd,lalim_conv,wght_th, &
2614	& ev, ep,epmlmMm,eplaMm, &
2615	& wdtrainA,wdtrainM)
2616	! RomP <<<
2617
2618	!IM begin
2619	!c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2620	!c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2621	!IM end
2622	!IM cf. FH
2623	clwcon0=qcondc
2624	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2625
2626	do i = 1, klon
2627	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2628	enddo
2629
2630	ELSE ! ok_cvl
2631
2632	!c MAF conema3 ne contient pas les traceurs
2633	CALL conema3 (dtime, &
2634	& paprs,pplay,t_seri,q_seri, &
2635	& u_seri,v_seri,tr_seri,ntra, &
2636	& ema_work1,ema_work2, &
2637	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2638	& rain_con, snow_con, ibas_con, itop_con, &
2639	& upwd,dnwd,dnwd0,bas,top, &
2640	& Ma,cape,tvp,rflag, &
2641	& pbase &
2642	& ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2643	& ,clwcon0)
2644
2645	ENDIF ! ok_cvl
2646
2647	!c
2648	!c Correction precip
2649	rain_con = rain_con * cvl_corr
2650	snow_con = snow_con * cvl_corr
2651	!c
2652
2653	IF (.NOT. ok_gust) THEN
2654	do i = 1, klon
2655	wd(i)=0.0
2656	enddo
2657	ENDIF
2658
2659	!c =================================================================== c
2660	!c Calcul des proprietes des nuages convectifs
2661	!c
2662
2663	!c calcul des proprietes des nuages convectifs
2664	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2665	call clouds_gno &
2666	& (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2667
2668	!c =================================================================== c
2669
2670	DO i = 1, klon
2671	itop_con(i) = min(max(itop_con(i),1),klev)
2672	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2673	ENDDO
2674
2675	DO i = 1, klon
2676	ema_pcb(i) = paprs(i,ibas_con(i))
2677	ENDDO
2678	DO i = 1, klon
2679	! L'idicage de itop_con peut cacher un pb potentiel
2680	! FH sous la dictee de JYG, CR
2681	ema_pct(i) = paprs(i,itop_con(i)+1)
2682
2683	if (itop_con(i).gt.klev-3) then
2684	if(prt_level >= 9) then
2685	write(lunout,*)'La convection monte trop haut '
2686	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2687	endif
2688	endif
2689	ENDDO
2690	ELSE IF (iflag_con.eq.0) THEN
2691	write(lunout,*) 'On n appelle pas la convection'
2692	clwcon0=0.
2693	rnebcon0=0.
2694	d_t_con=0.
2695	d_q_con=0.
2696	d_u_con=0.
2697	d_v_con=0.
2698	rain_con=0.
2699	snow_con=0.
2700	bas=1
2701	top=1
2702	ELSE
2703	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2704	CALL abort
2705	ENDIF
2706
2707	!c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2708	!c . d_u_con, d_v_con)
2709
2710	!-----------------------------------------------------------------------------------------
2711	! ajout des tendances de la diffusion turbulente
2712	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2713	!-----------------------------------------------------------------------------------------
2714
2715	if (mydebug) then
2716	call writefield_phy('u_seri',u_seri,llm)
2717	call writefield_phy('v_seri',v_seri,llm)
2718	call writefield_phy('t_seri',t_seri,llm)
2719	call writefield_phy('q_seri',q_seri,llm)
2720	endif
2721
2722	!IM
2723	IF (ip_ebil_phy.ge.2) THEN
2724	ztit='after convect'
2725	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2726	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2727	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2728	call diagphy(airephy,ztit,ip_ebil_phy &
2729	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2730	& , zero_v, rain_con, snow_con, ztsol &
2731	& , d_h_vcol, d_qt, d_ec &
2732	& , fs_bound, fq_bound )
2733	END IF
2734	!C
2735	PRINT *,' Lluis 4 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2736	fname = 'After convection'
2737	varnamechk = 'paprs'
2738	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
2739	varnamechk = 'pplay'
2740	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
2741	varnamechk = 'pphi'
2742	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
2743	varnamechk = 't_seri'
2744	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
2745	varnamechk = 'u_seri'
2746	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
2747	varnamechk = 'v_seri'
2748	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
2749	varnamechk = 'q_seri'
2750	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
2751
2752	IF (check) THEN
2753	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2754	WRITE(lunout,*)"aprescon=", za
2755	zx_t = 0.0
2756	za = 0.0
2757	DO i = 1, klon
2758	za = za + airephy(i)/REAL(klon)
2759	zx_t = zx_t + (rain_con(i)+ &
2760	& snow_con(i))*airephy(i)/REAL(klon)
2761	ENDDO
2762	zx_t = zx_t/za*dtime
2763	WRITE(lunout,*)"Precip=", zx_t
2764	ENDIF
2765	IF (zx_ajustq) THEN
2766	DO i = 1, klon
2767	z_apres(i) = 0.0
2768	ENDDO
2769	DO k = 1, klev
2770	DO i = 1, klon
2771	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2772	& *(paprs(i,k)-paprs(i,k+1))/RG
2773	ENDDO
2774	ENDDO
2775	DO i = 1, klon
2776	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2777	& /z_apres(i)
2778	ENDDO
2779	DO k = 1, klev
2780	DO i = 1, klon
2781	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2782	& z_factor(i).LT.(1.0-1.0E-08)) THEN
2783	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2784	ENDIF
2785	ENDDO
2786	ENDDO
2787	ENDIF
2788	zx_ajustq=.FALSE.
2789
2790	!c
2791	!c=============================================================================
2792	!cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2793	!cpour la couche limite diffuse pour l instant
2794	!c
2795	if (iflag_wake>=1) then
2796	DO k=1,klev
2797	DO i=1,klon
2798	dt_dwn(i,k) = ftd(i,k)
2799	wdt_PBL(i,k) = 0.
2800	dq_dwn(i,k) = fqd(i,k)
2801	wdq_PBL(i,k) = 0.
2802	M_dwn(i,k) = dnwd0(i,k)
2803	M_up(i,k) = upwd(i,k)
2804	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2805	udt_PBL(i,k) = 0.
2806	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2807	udq_PBL(i,k) = 0.
2808	ENDDO
2809	ENDDO
2810
2811	if (iflag_wake==2) then
2812	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2813	DO k = 1,klev
2814	dt_dwn(:,k)= dt_dwn(:,k)+ &
2815	& ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2816	dq_dwn(:,k)= dq_dwn(:,k)+ &
2817	& ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2818	ENDDO
2819	endif
2820	!c
2821	!ccalcul caracteristiques de la poche froide
2822	call calWAKE (paprs,pplay,dtime &
2823	& ,t_seri,q_seri,omega &
2824	& ,dt_dwn,dq_dwn,M_dwn,M_up &
2825	& ,dt_a,dq_a,sigd &
2826	& ,wdt_PBL,wdq_PBL &
2827	& ,udt_PBL,udq_PBL &
2828	& ,wake_deltat,wake_deltaq,wake_dth &
2829	& ,wake_h,wake_s,wake_dens &
2830	& ,wake_pe,wake_fip,wake_gfl &
2831	& ,dt_wake,dq_wake &
2832	& ,wake_k, t_undi,q_undi &
2833	& ,wake_omgbdth,wake_dp_omgb &
2834	& ,wake_dtKE,wake_dqKE &
2835	& ,wake_dtPBL,wake_dqPBL &
2836	& ,wake_omg,wake_dp_deltomg &
2837	& ,wake_spread,wake_Cstar,wake_d_deltat_gw &
2838	& ,wake_ddeltat,wake_ddeltaq)
2839	!c
2840	!-----------------------------------------------------------------------------------------
2841	! ajout des tendances des poches froides
2842	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2843	! coherent avec les autres d_t_...
2844	d_t_wake(:,:)=dt_wake(:,:)*dtime
2845	d_q_wake(:,:)=dq_wake(:,:)*dtime
2846	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2847	!-----------------------------------------------------------------------------------------
2848
2849	endif
2850	!c
2851	!c===================================================================
2852	!cJYG
2853	IF (ip_ebil_phy.ge.2) THEN
2854	ztit='after wake'
2855	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2856	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2857	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2858	call diagphy(airephy,ztit,ip_ebil_phy &
2859	& , zero_v, zero_v, zero_v, zero_v, zero_v &
2860	& , zero_v, zero_v, zero_v, ztsol &
2861	& , d_h_vcol, d_qt, d_ec &
2862	& , fs_bound, fq_bound )
2863	END IF
2864	PRINT *,' Lluis 5 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2865	fname = 'After wakes'
2866	varnamechk = 'paprs'
2867	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
2868	varnamechk = 'pplay'
2869	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
2870	varnamechk = 'pphi'
2871	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
2872	varnamechk = 't_seri'
2873	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
2874	varnamechk = 'u_seri'
2875	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
2876	varnamechk = 'v_seri'
2877	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
2878	varnamechk = 'q_seri'
2879	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
2880
2881	!c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2882	!c
2883	!c===================================================================
2884	!c Convection seche (thermiques ou ajustement)
2885	!c===================================================================
2886	!c
2887	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2888	& ,seuil_inversion,weak_inversion,dthmin)
2889
2890
2891
2892	d_t_ajsb(:,:)=0.
2893	d_q_ajsb(:,:)=0.
2894	d_t_ajs(:,:)=0.
2895	d_u_ajs(:,:)=0.
2896	d_v_ajs(:,:)=0.
2897	d_q_ajs(:,:)=0.
2898	clwcon0th(:,:)=0.
2899	!c
2900	!c fm_therm(:,:)=0.
2901	!c entr_therm(:,:)=0.
2902	!c detr_therm(:,:)=0.
2903	!c
2904
2905	IF(prt_level>9)WRITE(lunout,*) &
2906	& 'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2907	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2908	if(iflag_thermals.lt.0) then
2909	!c Rien
2910	!c ====
2911	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2912
2913
2914	else
2915
2916	!c Thermiques
2917	!c ==========
2918	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2919	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2920
2921
2922	!ccc nrlmd le 10/04/2012
2923	DO k=1,klev+1
2924	DO i=1,klon
2925	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2926	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2927	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2928	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2929	ENDDO
2930	ENDDO
2931	!ccc fin nrlmd le 10/04/2012
2932
2933	if (iflag_thermals>=1) then
2934	call calltherm(pdtphys &
2935	& ,pplay,paprs,pphi,weak_inversion &
2936	& ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
2937	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2938	& ,fm_therm,entr_therm,detr_therm &
2939	& ,zqasc,clwcon0th,lmax_th,ratqscth &
2940	& ,ratqsdiff,zqsatth &
2941	!con rajoute ale et alp, et les caracteristiques de la couche alim &
2942	& ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2943	& ,ztv,zpspsk,ztla,zthl &
2944	!ccc nrlmd le 10/04/2012 &
2945	& ,pbl_tke_input,pctsrf,omega,airephy &
2946	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2947	& ,n2,s2,ale_bl_stat &
2948	& ,therm_tke_max,env_tke_max &
2949	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2950	& ,alp_bl_conv,alp_bl_stat &
2951	!ccc fin nrlmd le 10/04/2012 &
2952	& ,zqla,ztva )
2953
2954	!ccc nrlmd le 10/04/2012
2955	!c-----------Stochastic triggering-----------
2956	if (iflag_trig_bl.ge.1) then
2957	!c
2958	IF (prt_level .GE. 10) THEN
2959	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2960	& cin, ale_bl_stat, alp_bl_stat
2961	ENDIF
2962
2963	!c----Initialisations
2964	do i=1,klon
2965	proba_notrig(i)=1.
2966	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2967	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2968	tau_trig(i)=tau_trig_shallow
2969	else
2970	tau_trig(i)=tau_trig_deep
2971	endif
2972	enddo
2973	!c
2974	IF (prt_level .GE. 10) THEN
2975	print *,'random_notrig, tau_trig ', &
2976	& random_notrig, tau_trig
2977	print *,'s_trig,s2,n2 ', &
2978	& s_trig,s2,n2
2979	ENDIF
2980
2981	!c----Tirage al\'eatoire et calcul de ale_bl_trig
2982	do i=1,klon
2983	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2984	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2985	& (n2(i)*dtime/tau_trig(i))
2986	!c print *, 'proba_notrig(i) ',proba_notrig(i)
2987	if (random_notrig(i) .ge. proba_notrig(i)) then
2988	ale_bl_trig(i)=ale_bl_stat(i)
2989	else
2990	ale_bl_trig(i)=0.
2991	endif
2992	else
2993	proba_notrig(i)=1.
2994	random_notrig(i)=0.
2995	ale_bl_trig(i)=0.
2996	endif
2997	enddo
2998	!c
2999	IF (prt_level .GE. 10) THEN
3000	print *,'proba_notrig, ale_bl_trig ', &
3001	& proba_notrig, ale_bl_trig
3002	ENDIF
3003
3004	endif !(iflag_trig_bl)
3005
3006	!c-----------Statistical closure-----------
3007	if (iflag_clos_bl.ge.1) then
3008
3009	do i=1,klon
3010	alp_bl(i)=alp_bl_stat(i)
3011	enddo
3012
3013	else
3014
3015	alp_bl_stat(:)=0.
3016
3017	endif !(iflag_clos_bl)
3018
3019	IF (prt_level .GE. 10) THEN
3020	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
3021	ENDIF
3022
3023	!ccc fin nrlmd le 10/04/2012
3024
3025	! ----------------------------------------------------------------------
3026	! Transport de la TKE par les panaches thermiques.
3027	! FH : 2010/02/01
3028	! if (iflag_pbl.eq.10) then
3029	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3030	! s rg,paprs,pbl_tke)
3031	! endif
3032	! ----------------------------------------------------------------------
3033	!IM/FH: 2011/02/23
3034	! Couplage Thermiques/Emanuel seulement si T<0
3035	if (iflag_coupl==2) then
3036	print*,'Couplage Thermiques/Emanuel seulement si T<0'
3037	do i=1,klon
3038	if (t_seri(i,lmax_th(i))>273.) then
3039	Ale_bl(i)=0.
3040	endif
3041	enddo
3042	endif
3043
3044	do i=1,klon
3045	zmax_th(i)=pphi(i,lmax_th(i))/rg
3046	enddo
3047
3048	endif
3049
3050
3051	!c Ajustement sec
3052	!c ==============
3053
3054	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3055	! a partir du sommet des thermiques.
3056	! Dans le cas contraire, on demarre au niveau 1.
3057
3058	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
3059
3060	if(iflag_thermals.eq.0) then
3061	IF(prt_level>9)WRITE(lunout,*)'ajsec'
3062	limbas(:)=1
3063	else
3064	limbas(:)=lmax_th(:)
3065	endif
3066
3067	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3068	! pour des test de convergence numerique.
3069	! Le nouveau ajsec est a priori mieux, meme pour le cas
3070	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3071	! non nulles numeriquement pour des mailles non concernees.
3072
3073	if (iflag_thermals.eq.0) then
3074	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3075	& , d_t_ajsb, d_q_ajsb)
3076	else
3077	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3078	& , d_t_ajsb, d_q_ajsb)
3079	endif
3080
3081	!-----------------------------------------------------------------------------------------
3082	! ajout des tendances de l'ajustement sec ou des thermiques
3083	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
3084	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3085	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3086
3087	!-----------------------------------------------------------------------------------------
3088
3089	endif
3090
3091	endif
3092
3093	!c
3094	!c===================================================================
3095	!IM
3096	IF (ip_ebil_phy.ge.2) THEN
3097	ztit='after dry_adjust'
3098	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3099	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3100	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3101	call diagphy(airephy,ztit,ip_ebil_phy &
3102	& , zero_v, zero_v, zero_v, zero_v, zero_v &
3103	& , zero_v, zero_v, zero_v, ztsol &
3104	& , d_h_vcol, d_qt, d_ec &
3105	& , fs_bound, fq_bound )
3106	END IF
3107	PRINT *,' Lluis 6 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3108	fname = 'After thermals'
3109	varnamechk = 'paprs'
3110	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
3111	varnamechk = 'pplay'
3112	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
3113	varnamechk = 'pphi'
3114	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
3115	varnamechk = 't_seri'
3116	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
3117	varnamechk = 'u_seri'
3118	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
3119	varnamechk = 'v_seri'
3120	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
3121	varnamechk = 'q_seri'
3122	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
3123
3124	!c-------------------------------------------------------------------------
3125	! Computation of ratqs, the width (normalized) of the subrid scale
3126	! water distribution
3127	CALL calcratqs(klon,klev,prt_level,lunout, &
3128	& iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
3129	& ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
3130	& ptconv,ptconvth,clwcon0th, rnebcon0th, &
3131	& paprs,pplay,q_seri,zqsat,fm_therm, &
3132	& ratqs,ratqsc)
3133
3134
3135	!c
3136	!c Appeler le processus de condensation a grande echelle
3137	!c et le processus de precipitation
3138	!c-------------------------------------------------------------------------
3139	IF (prt_level .GE.10) THEN
3140	print *,' ->fisrtilp '
3141	ENDIF
3142	!c-------------------------------------------------------------------------
3143	CALL fisrtilp(dtime,paprs,pplay, &
3144	& t_seri, q_seri,ptconv,ratqs, &
3145	& d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
3146	& rain_lsc, snow_lsc, &
3147	& pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3148	& frac_impa, frac_nucl, beta_prec_fisrt, &
3149	& prfl, psfl, rhcl, &
3150	& zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
3151
3152	WHERE (rain_lsc < 0) rain_lsc = 0.
3153	WHERE (snow_lsc < 0) snow_lsc = 0.
3154	!-----------------------------------------------------------------------------------------
3155	! ajout des tendances de la diffusion turbulente
3156	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
3157	!-----------------------------------------------------------------------------------------
3158	DO k = 1, klev
3159	DO i = 1, klon
3160	cldfra(i,k) = rneb(i,k)
3161	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3162	ENDDO
3163	ENDDO
3164	IF (check) THEN
3165	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
3166	WRITE(lunout,*)"apresilp=", za
3167	zx_t = 0.0
3168	za = 0.0
3169	DO i = 1, klon
3170	za = za + airephy(i)/REAL(klon)
3171	zx_t = zx_t + (rain_lsc(i) &
3172	& + snow_lsc(i))*airephy(i)/REAL(klon)
3173	ENDDO
3174	zx_t = zx_t/za*dtime
3175	WRITE(lunout,*)"Precip=", zx_t
3176	ENDIF
3177	!IM
3178	IF (ip_ebil_phy.ge.2) THEN
3179	ztit='after fisrt'
3180	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3181	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3182	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3183	call diagphy(airephy,ztit,ip_ebil_phy &
3184	& , zero_v, zero_v, zero_v, zero_v, zero_v &
3185	& , zero_v, rain_lsc, snow_lsc, ztsol &
3186	& , d_h_vcol, d_qt, d_ec &
3187	& , fs_bound, fq_bound )
3188	END IF
3189	PRINT *,' Lluis 7 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3190	fname = 'After firstilp'
3191	varnamechk = 'paprs'
3192	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
3193	varnamechk = 'pplay'
3194	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
3195	varnamechk = 'pphi'
3196	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
3197	varnamechk = 't_seri'
3198	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
3199	varnamechk = 'u_seri'
3200	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
3201	varnamechk = 'v_seri'
3202	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
3203	varnamechk = 'q_seri'
3204	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
3205
3206	if (mydebug) then
3207	call writefield_phy('u_seri',u_seri,llm)
3208	call writefield_phy('v_seri',v_seri,llm)
3209	call writefield_phy('t_seri',t_seri,llm)
3210	call writefield_phy('q_seri',q_seri,llm)
3211	endif
3212	!c
3213	!c-------------------------------------------------------------------
3214	!c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3215	!c-------------------------------------------------------------------
3216
3217	!c 1. NUAGES CONVECTIFS
3218	!c
3219	!IM cf FH
3220	!c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
3221	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
3222	snow_tiedtke=0.
3223	!c print*,'avant calcul de la pseudo precip '
3224	!c print*,'iflag_cldcon',iflag_cldcon
3225	if (iflag_cldcon.eq.-1) then
3226	rain_tiedtke=rain_con
3227	else
3228	!c print*,'calcul de la pseudo precip '
3229	rain_tiedtke=0.
3230	!c print*,'calcul de la pseudo precip 0'
3231	do k=1,klev
3232	do i=1,klon
3233	if (d_q_con(i,k).lt.0.) then
3234	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3235	& *(paprs(i,k)-paprs(i,k+1))/rg
3236	endif
3237	enddo
3238	enddo
3239	endif
3240	!c
3241	!c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3242	!c
3243
3244	!c Nuages diagnostiques pour Tiedtke
3245	CALL diagcld1(paprs,pplay, &
3246	!IM cf FH . rain_con,snow_con,ibas_con,itop_con, &
3247	& rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3248	& diafra,dialiq)
3249	DO k = 1, klev
3250	DO i = 1, klon
3251	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3252	cldliq(i,k) = dialiq(i,k)
3253	cldfra(i,k) = diafra(i,k)
3254	ENDIF
3255	ENDDO
3256	ENDDO
3257
3258	ELSE IF (iflag_cldcon.ge.3) THEN
3259	!c On prend pour les nuages convectifs le max du calcul de la
3260	!c convection et du calcul du pas de temps precedent diminue d'un facteur
3261	!c facttemps
3262	facteur = pdtphys *facttemps
3263	do k=1,klev
3264	do i=1,klon
3265	rnebcon(i,k)=rnebcon(i,k)*facteur
3266	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
3267	& then
3268	rnebcon(i,k)=rnebcon0(i,k)
3269	clwcon(i,k)=clwcon0(i,k)
3270	endif
3271	enddo
3272	enddo
3273
3274	!c
3275	!cjq - introduce the aerosol direct and first indirect radiative forcings
3276	!cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3277	IF (flag_aerosol .gt. 0) THEN
3278	IF (.NOT. aerosol_couple) &
3279	& CALL readaerosol_optic( &
3280	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3281	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3282	& mass_solu_aero, mass_solu_aero_pi, &
3283	& tau_aero, piz_aero, cg_aero, &
3284	& tausum_aero, tau3d_aero)
3285	ELSE
3286	tausum_aero(:,:,:) = 0.
3287	tau_aero(:,:,:,:) = 0.
3288	piz_aero(:,:,:,:) = 0.
3289	cg_aero(:,:,:,:) = 0.
3290	ENDIF
3291	!c
3292	!c--STRAT AEROSOL
3293	!c--updates tausum_aero,tau_aero,piz_aero,cg_aero
3294	IF (flag_aerosol_strat) THEN
3295	PRINT *,'appel a readaerosolstrat', mth_cur
3296	CALL readaerosolstrato(debut)
3297	ENDIF
3298	!c--fin STRAT AEROSOL
3299
3300	!IM calcul nuages par le simulateur ISCCP
3301	!c
3302	#ifdef histISCCP
3303	IF (ok_isccp) THEN
3304	!c
3305	!IM lecture invtau, tautab des fichiers formattes
3306	!c
3307	IF (debut) THEN
3308	!$OMP MASTER
3309	!c
3310	open(99,file='tautab.formatted', FORM='FORMATTED')
3311	read(99,'(f30.20)') tautab_omp
3312	close(99)
3313	!c
3314	open(99,file='invtau.formatted',form='FORMATTED')
3315	read(99,'(i10)') invtau_omp
3316
3317	!c print*,'calcul_simulISCCP invtau_omp',invtau_omp
3318	!c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
3319
3320	close(99)
3321	!$OMP END MASTER
3322	!$OMP BARRIER
3323	tautab=tautab_omp
3324	invtau=invtau_omp
3325	!c
3326	ENDIF !debut
3327	!c
3328	!IM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
3329	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
3330	#include "calcul_simulISCCP.h"
3331	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
3332	ENDIF !ok_isccp
3333	#endif
3334
3335	!c On prend la somme des fractions nuageuses et des contenus en eau
3336
3337	if (iflag_cldcon>=5) then
3338
3339	do k=1,klev
3340	ptconvth(:,k)=fm_therm(:,k+1)>0.
3341	enddo
3342
3343	if (iflag_coupl==4) then
3344
3345	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3346	! convectives et lsc dans la partie des thermiques
3347	! Le controle par iflag_coupl est peut etre provisoire.
3348	do k=1,klev
3349	do i=1,klon
3350	if (ptconv(i,k).and.ptconvth(i,k)) then
3351	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3352	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3353	else if (ptconv(i,k)) then
3354	cldfra(i,k)=rnebcon(i,k)
3355	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3356	endif
3357	enddo
3358	enddo
3359
3360	else if (iflag_coupl==5) then
3361	do k=1,klev
3362	do i=1,klon
3363	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3364	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3365	enddo
3366	enddo
3367
3368	else
3369
3370	! Si on est sur un point touche par la convection profonde et pas
3371	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3372	! de la convection profonde.
3373
3374	!IM/FH: 2011/02/23
3375	! definition des points sur lesquels ls thermiques sont actifs
3376
3377	do k=1,klev
3378	do i=1,klon
3379	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3380	cldfra(i,k)=rnebcon(i,k)
3381	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3382	endif
3383	enddo
3384	enddo
3385
3386	endif
3387
3388	else
3389
3390	! Ancienne version
3391	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3392	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3393	endif
3394
3395	ENDIF
3396
3397	! plulsc(:)=0.
3398	! do k=1,klev,-1
3399	! do i=1,klon
3400	! zzz=prfl(:,k)+psfl(:,k)
3401	! if (.not.ptconvth.zzz.gt.0.)
3402	! enddo prfl, psfl,
3403	! enddo
3404	!c
3405	!c 2. NUAGES STARTIFORMES
3406	!c
3407	IF (ok_stratus) THEN
3408	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3409	DO k = 1, klev
3410	DO i = 1, klon
3411	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3412	cldliq(i,k) = dialiq(i,k)
3413	cldfra(i,k) = diafra(i,k)
3414	ENDIF
3415	ENDDO
3416	ENDDO
3417	ENDIF
3418	!c
3419	!c Precipitation totale
3420	!c
3421	DO i = 1, klon
3422	rain_fall(i) = rain_con(i) + rain_lsc(i)
3423	snow_fall(i) = snow_con(i) + snow_lsc(i)
3424	ENDDO
3425	!IM
3426	IF (ip_ebil_phy.ge.2) THEN
3427	ztit="after diagcld"
3428	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3429	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3430	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3431	call diagphy(airephy,ztit,ip_ebil_phy &
3432	& , zero_v, zero_v, zero_v, zero_v, zero_v &
3433	& , zero_v, zero_v, zero_v, ztsol &
3434	& , d_h_vcol, d_qt, d_ec &
3435	& , fs_bound, fq_bound )
3436	END IF
3437	PRINT *,' Lluis 8 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3438	fname = 'After coupling convection+thermics+wakes'
3439	varnamechk = 'paprs'
3440	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
3441	varnamechk = 'pplay'
3442	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
3443	varnamechk = 'pphi'
3444	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
3445	varnamechk = 't_seri'
3446	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
3447	varnamechk = 'u_seri'
3448	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
3449	varnamechk = 'v_seri'
3450	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
3451	varnamechk = 'q_seri'
3452	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
3453
3454	!c
3455	!c Calculer l'humidite relative pour diagnostique
3456	!c
3457	DO k = 1, klev
3458	DO i = 1, klon
3459	zx_t = t_seri(i,k)
3460	IF (thermcep) THEN
3461	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3462	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3463	zx_qs = MIN(0.5,zx_qs)
3464	zcor = 1./(1.-retv*zx_qs)
3465	zx_qs = zx_qs*zcor
3466	ELSE
3467	IF (zx_t.LT.t_coup) THEN
3468	zx_qs = qsats(zx_t)/pplay(i,k)
3469	ELSE
3470	zx_qs = qsatl(zx_t)/pplay(i,k)
3471	ENDIF
3472	ENDIF
3473	zx_rh(i,k) = q_seri(i,k)/zx_qs
3474	zqsat(i,k)=zx_qs
3475	ENDDO
3476	ENDDO
3477
3478	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3479	!c equivalente a 2m (tpote) pour diagnostique
3480	!c
3481	DO i = 1, klon
3482	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3483	IF (thermcep) THEN
3484	IF(zt2m(i).LT.RTT) then
3485	Lheat=RLSTT
3486	ELSE
3487	Lheat=RLVTT
3488	ENDIF
3489	ELSE
3490	IF (zt2m(i).LT.RTT) THEN
3491	Lheat=RLSTT
3492	ELSE
3493	Lheat=RLVTT
3494	ENDIF
3495	ENDIF
3496	tpote(i) = tpot(i)* &
3497	& EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3498	ENDDO
3499
3500	IF (type_trac == 'inca') THEN
3501	#ifdef INCA
3502	CALL VTe(VTphysiq)
3503	CALL VTb(VTinca)
3504	calday = REAL(days_elapsed + 1) + jH_cur
3505
3506	call chemtime(itap+itau_phy-1, date0, dtime)
3507	IF (config_inca == 'aero') THEN
3508	CALL AEROSOL_METEO_CALC( &
3509	& calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3510	& prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3511	END IF
3512
3513	zxsnow_dummy(:) = 0.0
3514
3515	CALL chemhook_begin (calday, &
3516	& days_elapsed+1, &
3517	& jH_cur, &
3518	& pctsrf(1,1), &
3519	& rlat, &
3520	& rlon, &
3521	& airephy, &
3522	& paprs, &
3523	& pplay, &
3524	& coefh(:,:,is_ave), &
3525	& pphi, &
3526	& t_seri, &
3527	& u, &
3528	& v, &
3529	& wo(:, :, 1), &
3530	& q_seri, &
3531	& zxtsol, &
3532	& zxsnow_dummy, &
3533	& solsw, &
3534	& albsol1, &
3535	& rain_fall, &
3536	& snow_fall, &
3537	& itop_con, &
3538	& ibas_con, &
3539	& cldfra, &
3540	& iim, &
3541	& jjm, &
3542	& tr_seri, &
3543	& ftsol, &
3544	& paprs, &
3545	& cdragh, &
3546	& cdragm, &
3547	& pctsrf, &
3548	& pdtphys, &
3549	& itap)
3550
3551	CALL VTe(VTinca)
3552	CALL VTb(VTphysiq)
3553	#endif
3554	END IF !type_trac = inca
3555	!c
3556	!c Calculer les parametres optiques des nuages et quelques
3557	!c parametres pour diagnostiques:
3558	!c
3559
3560	IF (aerosol_couple) THEN
3561	mass_solu_aero(:,:) = ccm(:,:,1)
3562	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3563	END IF
3564
3565	if (ok_newmicro) then
3566	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3567	& paprs, pplay, t_seri, cldliq, cldfra, &
3568	& cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3569	& flwp, fiwp, flwc, fiwc, &
3570	& mass_solu_aero, mass_solu_aero_pi, &
3571	& cldtaupi, re, fl, ref_liq, ref_ice)
3572	else
3573	CALL nuage (paprs, pplay, &
3574	& t_seri, cldliq, cldfra, cldtau, cldemi, &
3575	& cldh, cldl, cldm, cldt, cldq, &
3576	& ok_aie, &
3577	& mass_solu_aero, mass_solu_aero_pi, &
3578	& bl95_b0, bl95_b1, &
3579	& cldtaupi, re, fl)
3580	endif
3581	!c
3582	!IM betaCRF
3583	!c
3584	cldtaurad = cldtau
3585	cldtaupirad = cldtaupi
3586	cldemirad = cldemi &
3587	&!c
3588	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3589	&lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3590	!c
3591	!c global
3592	!c
3593	DO k=1, klev
3594	DO i=1, klon
3595	if (pplay(i,k).GE.pfree) THEN
3596	beta(i,k) = beta_pbl
3597	else
3598	beta(i,k) = beta_free
3599	endif
3600	if (mskocean_beta) THEN
3601	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3602	endif
3603	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3604	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3605	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3606	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3607	ENDDO
3608	ENDDO
3609	!c
3610	else
3611	!c
3612	!c regional
3613	!c
3614	DO k=1, klev
3615	DO i=1,klon
3616	!c
3617	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3618	& rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3619	if (pplay(i,k).GE.pfree) THEN
3620	beta(i,k) = beta_pbl
3621	else
3622	beta(i,k) = beta_free
3623	endif
3624	if (mskocean_beta) THEN
3625	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3626	endif
3627	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3628	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3629	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3630	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3631	endif
3632	!c
3633	ENDDO
3634	ENDDO
3635	!c
3636	endif
3637	!c
3638	!c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3639	!c
3640	IF (MOD(itaprad,radpas).EQ.0) THEN
3641
3642	DO i = 1, klon
3643	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3644	& + falb1(i,is_lic) * pctsrf(i,is_lic) &
3645	& + falb1(i,is_ter) * pctsrf(i,is_ter) &
3646	& + falb1(i,is_sic) * pctsrf(i,is_sic)
3647	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3648	& + falb2(i,is_lic) * pctsrf(i,is_lic) &
3649	& + falb2(i,is_ter) * pctsrf(i,is_ter) &
3650	& + falb2(i,is_sic) * pctsrf(i,is_sic)
3651	ENDDO
3652
3653	if (mydebug) then
3654	call writefield_phy('u_seri',u_seri,llm)
3655	call writefield_phy('v_seri',v_seri,llm)
3656	call writefield_phy('t_seri',t_seri,llm)
3657	call writefield_phy('q_seri',q_seri,llm)
3658	endif
3659
3660	IF (aerosol_couple) THEN
3661	#ifdef INCA
3662	CALL radlwsw_inca &
3663	& (kdlon,kflev,dist, rmu0, fract, solaire, &
3664	& paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3665	& wo(:, :, 1), &
3666	& cldfrarad, cldemirad, cldtaurad, &
3667	& heat,heat0,cool,cool0,radsol,albpla, &
3668	& topsw,toplw,solsw,sollw, &
3669	& sollwdown, &
3670	& topsw0,toplw0,solsw0,sollw0, &
3671	& lwdn0, lwdn, lwup0, lwup, &
3672	& swdn0, swdn, swup0, swup, &
3673	& ok_ade, ok_aie, &
3674	& tau_aero, piz_aero, cg_aero, &
3675	& topswad_aero, solswad_aero, &
3676	& topswad0_aero, solswad0_aero, &
3677	& topsw_aero, topsw0_aero, &
3678	& solsw_aero, solsw0_aero, &
3679	& cldtaupirad, &
3680	& topswai_aero, solswai_aero)
3681
3682	#endif
3683	ELSE
3684	!c
3685	!IM calcul radiatif pour le cas actuel
3686	!c
3687	RCO2 = RCO2_act
3688	RCH4 = RCH4_act
3689	RN2O = RN2O_act
3690	RCFC11 = RCFC11_act
3691	RCFC12 = RCFC12_act
3692	!c
3693	IF (prt_level .GE.10) THEN
3694	print *,' ->radlwsw, number 1 '
3695	ENDIF
3696	!c
3697	CALL radlwsw &
3698	& (dist, rmu0, fract, &
3699	& paprs, pplay,zxtsol,albsol1, albsol2, &
3700	& t_seri,q_seri,wo, &
3701	& cldfrarad, cldemirad, cldtaurad, &
3702	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3703	& flag_aerosol_strat, &
3704	& tau_aero, piz_aero, cg_aero, &
3705	& cldtaupirad,new_aod, &
3706	& zqsat, flwc, fiwc, &
3707	& heat,heat0,cool,cool0,radsol,albpla, &
3708	& topsw,toplw,solsw,sollw, &
3709	& sollwdown, &
3710	& topsw0,toplw0,solsw0,sollw0, &
3711	& lwdn0, lwdn, lwup0, lwup, &
3712	& swdn0, swdn, swup0, swup, &
3713	& topswad_aero, solswad_aero, &
3714	& topswai_aero, solswai_aero, &
3715	& topswad0_aero, solswad0_aero, &
3716	& topsw_aero, topsw0_aero, &
3717	& solsw_aero, solsw0_aero, &
3718	& topswcf_aero, solswcf_aero)
3719
3720	!c
3721	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3722	!IM des taux doit etre different du taux actuel
3723	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3724	!c
3725
3726	if (ok_4xCO2atm) then
3727	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3728	&RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3729	&RCFC12_per.NE.RCFC12_act) THEN
3730	!c
3731	RCO2 = RCO2_per
3732	RCH4 = RCH4_per
3733	RN2O = RN2O_per
3734	RCFC11 = RCFC11_per
3735	RCFC12 = RCFC12_per
3736	!c
3737	IF (prt_level .GE.10) THEN
3738	print *,' ->radlwsw, number 2 '
3739	ENDIF
3740	!c
3741	CALL radlwsw &
3742	& (dist, rmu0, fract, &
3743	& paprs, pplay,zxtsol,albsol1, albsol2, &
3744	& t_seri,q_seri,wo, &
3745	& cldfra, cldemi, cldtau, &
3746	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3747	& flag_aerosol_strat, &
3748	& tau_aero, piz_aero, cg_aero, &
3749	& cldtaupi,new_aod, &
3750	& zqsat, flwc, fiwc, &
3751	& heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3752	& topswp,toplwp,solswp,sollwp, &
3753	& sollwdownp, &
3754	& topsw0p,toplw0p,solsw0p,sollw0p, &
3755	& lwdn0p, lwdnp, lwup0p, lwupp, &
3756	& swdn0p, swdnp, swup0p, swupp, &
3757	& topswad_aerop, solswad_aerop, &
3758	& topswai_aerop, solswai_aerop, &
3759	& topswad0_aerop, solswad0_aerop, &
3760	& topsw_aerop, topsw0_aerop, &
3761	& solsw_aerop, solsw0_aerop, &
3762	& topswcf_aerop, solswcf_aerop)
3763	endif
3764	endif
3765	!c
3766	ENDIF ! aerosol_couple
3767	itaprad = 0
3768	ENDIF ! MOD(itaprad,radpas)
3769	itaprad = itaprad + 1
3770
3771	IF (iflag_radia.eq.0) THEN
3772	IF (prt_level.ge.9) THEN
3773	PRINT *,'--------------------------------------------------'
3774	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3775	PRINT *,'>>>> heat et cool mis a zero '
3776	PRINT *,'--------------------------------------------------'
3777	END IF
3778	heat=0.
3779	cool=0.
3780	sollw=0. ! MPL 01032011
3781	solsw=0.
3782	radsol=0.
3783	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3784	swup0=0.
3785	swdn=0.
3786	swdn0=0.
3787	lwup=0.
3788	lwup0=0.
3789	lwdn=0.
3790	lwdn0=0.
3791	END IF
3792
3793	!c
3794	!c Ajouter la tendance des rayonnements (tous les pas)
3795	!c
3796	DO k = 1, klev
3797	DO i = 1, klon
3798	t_seri(i,k) = t_seri(i,k) &
3799	& + (heat(i,k)-cool(i,k)) * dtime/RDAY
3800	ENDDO
3801	ENDDO
3802	!c
3803	if (mydebug) then
3804	call writefield_phy('u_seri',u_seri,llm)
3805	call writefield_phy('v_seri',v_seri,llm)
3806	call writefield_phy('t_seri',t_seri,llm)
3807	call writefield_phy('q_seri',q_seri,llm)
3808	endif
3809
3810	!IM
3811	IF (ip_ebil_phy.ge.2) THEN
3812	ztit='after rad'
3813	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3814	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3815	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3816	call diagphy(airephy,ztit,ip_ebil_phy &
3817	& , topsw, toplw, solsw, sollw, zero_v &
3818	& , zero_v, zero_v, zero_v, ztsol &
3819	& , d_h_vcol, d_qt, d_ec &
3820	& , fs_bound, fq_bound )
3821	END IF
3822	PRINT *,' Lluis 9 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3823	fname = 'After radiation'
3824	varnamechk = 'paprs'
3825	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
3826	varnamechk = 'pplay'
3827	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
3828	varnamechk = 'pphi'
3829	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
3830	varnamechk = 't_seri'
3831	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
3832	varnamechk = 'u_seri'
3833	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
3834	varnamechk = 'v_seri'
3835	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
3836	varnamechk = 'q_seri'
3837	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
3838
3839	!c
3840	!c
3841	!c Calculer l'hydrologie de la surface
3842	!c
3843	!c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3844	!c . agesno, ftsol,fqsurf,fsnow, ruis)
3845	!c
3846
3847	!c
3848	!c Calculer le bilan du sol et la derive de temperature (couplage)
3849	!c
3850	DO i = 1, klon
3851	!c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3852	!c a la demande de JLD
3853	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3854	ENDDO
3855	!c
3856	!cmoddeblott(jan95)
3857	!c Appeler le programme de parametrisation de l'orographie
3858	!c a l'echelle sous-maille:
3859	!c
3860	IF (prt_level .GE.10) THEN
3861	print *,' call orography ? ', ok_orodr
3862	ENDIF
3863	!c
3864	IF (ok_orodr) THEN
3865	!c
3866	!c selection des points pour lesquels le shema est actif:
3867	igwd=0
3868	DO i=1,klon
3869	itest(i)=0
3870	!c IF ((zstd(i).gt.10.0)) THEN
3871	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3872	itest(i)=1
3873	igwd=igwd+1
3874	idx(igwd)=i
3875	ENDIF
3876	ENDDO
3877	!c igwdim=MAX(1,igwd)
3878	!c
3879	IF (ok_strato) THEN
3880
3881	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3882	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3883	& igwd,idx,itest, &
3884	& t_seri, u_seri, v_seri, &
3885	& zulow, zvlow, zustrdr, zvstrdr, &
3886	& d_t_oro, d_u_oro, d_v_oro)
3887
3888	ELSE
3889	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3890	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3891	& igwd,idx,itest, &
3892	& t_seri, u_seri, v_seri, &
3893	& zulow, zvlow, zustrdr, zvstrdr, &
3894	& d_t_oro, d_u_oro, d_v_oro)
3895	ENDIF
3896	!c
3897	!c ajout des tendances
3898	!-----------------------------------------------------------------------------------------
3899	! ajout des tendances de la trainee de l'orographie
3900	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3901	!-----------------------------------------------------------------------------------------
3902	!c
3903	ENDIF ! fin de test sur ok_orodr
3904	!c
3905	if (mydebug) then
3906	call writefield_phy('u_seri',u_seri,llm)
3907	call writefield_phy('v_seri',v_seri,llm)
3908	call writefield_phy('t_seri',t_seri,llm)
3909	call writefield_phy('q_seri',q_seri,llm)
3910	endif
3911
3912	IF (ok_orolf) THEN
3913	!c
3914	!c selection des points pour lesquels le shema est actif:
3915	igwd=0
3916	DO i=1,klon
3917	itest(i)=0
3918	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3919	itest(i)=1
3920	igwd=igwd+1
3921	idx(igwd)=i
3922	ENDIF
3923	ENDDO
3924	!c igwdim=MAX(1,igwd)
3925	!c
3926	IF (ok_strato) THEN
3927
3928	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3929	& rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3930	& igwd,idx,itest, &
3931	& t_seri, u_seri, v_seri, &
3932	& zulow, zvlow, zustrli, zvstrli, &
3933	& d_t_lif, d_u_lif, d_v_lif )
3934
3935	ELSE
3936	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3937	& rlat,zmea,zstd,zpic, &
3938	& itest, &
3939	& t_seri, u_seri, v_seri, &
3940	& zulow, zvlow, zustrli, zvstrli, &
3941	& d_t_lif, d_u_lif, d_v_lif)
3942	ENDIF
3943	!c
3944	!-----------------------------------------------------------------------------------------
3945	! ajout des tendances de la portance de l'orographie
3946	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3947	!-----------------------------------------------------------------------------------------
3948	!c
3949	ENDIF ! fin de test sur ok_orolf
3950	!C HINES GWD PARAMETRIZATION
3951
3952	IF (ok_hines) then
3953
3954	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3955	& rlat,t_seri,u_seri,v_seri, &
3956	& zustrhi,zvstrhi, &
3957	& d_t_hin, d_u_hin, d_v_hin)
3958	!c
3959	!c ajout des tendances
3960	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
3961
3962	ENDIF
3963	!c
3964
3965	!c
3966	!IM cf. FLott BEG
3967	!C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3968
3969	if (mydebug) then
3970	call writefield_phy('u_seri',u_seri,llm)
3971	call writefield_phy('v_seri',v_seri,llm)
3972	call writefield_phy('t_seri',t_seri,llm)
3973	call writefield_phy('q_seri',q_seri,llm)
3974	endif
3975
3976	DO i = 1, klon
3977	zustrph(i)=0.
3978	zvstrph(i)=0.
3979	ENDDO
3980	DO k = 1, klev
3981	DO i = 1, klon
3982	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3983	& (paprs(i,k)-paprs(i,k+1))/rg
3984	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3985	& (paprs(i,k)-paprs(i,k+1))/rg
3986	ENDDO
3987	ENDDO
3988	!c
3989	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3990	!c
3991	! IF (is_sequential .and. ok_orodr) THEN
3992	! CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3993	! & ra,rg,romega, &
3994	! & rlat,rlon,pphis, &
3995	! & zustrdr,zustrli,zustrph, &
3996	! & zvstrdr,zvstrli,zvstrph, &
3997	! & paprs,u,v, &
3998	! & aam, torsfc)
3999	! ENDIF
4000	!IM cf. FLott END
4001	!IM
4002	IF (ip_ebil_phy.ge.2) THEN
4003	ztit='after orography'
4004	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
4005	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
4006	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
4007	call diagphy(airephy,ztit,ip_ebil_phy &
4008	& , zero_v, zero_v, zero_v, zero_v, zero_v &
4009	& , zero_v, zero_v, zero_v, ztsol &
4010	& , d_h_vcol, d_qt, d_ec &
4011	& , fs_bound, fq_bound )
4012	END IF
4013	PRINT *,' Lluis 10 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
4014	fname = 'after orography'
4015	varnamechk = 'paprs'
4016	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
4017	varnamechk = 'pplay'
4018	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
4019	varnamechk = 'pphi'
4020	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
4021	varnamechk = 't_seri'
4022	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
4023	varnamechk = 'u_seri'
4024	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
4025	varnamechk = 'v_seri'
4026	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
4027	varnamechk = 'q_seri'
4028	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
4029	!c
4030	!c
4031	!====================================================================
4032	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4033	!====================================================================
4034	! Abderrahmane 24.08.09
4035
4036	IF (ok_cosp) THEN
4037	! adeclarer
4038	#ifdef CPP_COSP
4039	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4040
4041	print*,'freq_cosp',freq_cosp
4042	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4043	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4044	! s ref_liq,ref_ice
4045	call phys_cosp(itap,dtime,freq_cosp, &
4046	& ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4047	& ecrit_mth,ecrit_day,ecrit_hf, &
4048	& klon,klev,rlon,rlat,presnivs,overlap, &
4049	& ref_liq,ref_ice, &
4050	& pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4051	& zu10m,zv10m,pphis, &
4052	& zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4053	& qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4054	& prfl(:,1:klev),psfl(:,1:klev), &
4055	& pmflxr(:,1:klev),pmflxs(:,1:klev), &
4056	& mr_ozone,cldtau, cldemi)
4057
4058	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4059	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4060	! M clMISR,
4061	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4062	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4063
4064	ENDIF
4065
4066	#endif
4067	ENDIF !ok_cosp
4068	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4069	!cAA
4070	!cAA Installation de l'interface online-offline pour traceurs
4071	!cAA
4072	!c====================================================================
4073	!c Calcul des tendances traceurs
4074	!c====================================================================
4075	!C
4076
4077	IF (type_trac=='repr') THEN
4078	sh_in(:,:) = q_seri(:,:)
4079	ELSE
4080	sh_in(:,:) = qx(:,:,ivap)
4081	END IF
4082
4083	call phytrac ( &
4084	& itap, days_elapsed+1, jH_cur, debut, &
4085	& lafin, dtime, u, v, t, &
4086	& paprs, pplay, pmfu, pmfd, &
4087	& pen_u, pde_u, pen_d, pde_d, &
4088	& cdragh, coefh(:,:,is_ave), fm_therm, entr_therm, &
4089	& u1, v1, ftsol, pctsrf, &
4090	& ustar, u10m, v10m, &
4091	& rlat, rlon, &
4092	& frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4093	& presnivs, pphis, pphi, albsol1, &
4094	& sh_in, rhcl, cldfra, rneb, &
4095	& diafra, cldliq, itop_con, ibas_con, &
4096	& pmflxr, pmflxs, prfl, psfl, &
4097	& da, phi, mp, upwd, &
4098	& phi2, d1a, dam, sij, &
4099	& wdtrainA, wdtrainM, sigd, clw,elij, &
4100	& ev, ep, epmlmMm, eplaMm, &
4101	& dnwd, aerosol_couple, flxmass_w, &
4102	& tau_aero, piz_aero, cg_aero, ccm, &
4103	& rfname, &
4104	& d_tr_dyn, &
4105	& tr_seri)
4106
4107	IF (offline) THEN
4108
4109	IF (prt_level.ge.9) &
4110	& print*,'Attention on met a 0 les thermiques pour phystoke'
4111	call phystokenc ( &
4112	& nlon,klev,pdtphys,rlon,rlat, &
4113	& t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4114	& fm_therm,entr_therm, &
4115	& cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf, &
4116	& frac_impa, frac_nucl, &
4117	& pphis,airephy,dtime,itap, &
4118	& qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4119
4120
4121	ENDIF
4122
4123	!c
4124	!c Calculer le transport de l'eau et de l'energie (diagnostique)
4125	!c
4126	CALL transp (paprs,zxtsol, &
4127	& t_seri, q_seri, u_seri, v_seri, zphi, &
4128	& ve, vq, ue, uq)
4129	!c
4130	!IM global posePB BEG
4131	IF(1.EQ.0) THEN
4132	!c
4133	CALL transp_lay (paprs,zxtsol, &
4134	& t_seri, q_seri, u_seri, v_seri, zphi, &
4135	& ve_lay, vq_lay, ue_lay, uq_lay)
4136	!c
4137	ENDIF !(1.EQ.0) THEN
4138	!IM global posePB END
4139	!c Accumuler les variables a stocker dans les fichiers histoire:
4140	!c
4141	fname = 'after phytrac'
4142	varnamechk = 'paprs'
4143	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
4144	varnamechk = 'pplay'
4145	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
4146	varnamechk = 'pphi'
4147	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
4148	varnamechk = 't_seri'
4149	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
4150	varnamechk = 'u_seri'
4151	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
4152	varnamechk = 'v_seri'
4153	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
4154	varnamechk = 'q_seri'
4155	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
4156
4157	!================================================================
4158	! Conversion of kinetic and potential energy into heat, for
4159	! parameterisation of subgrid-scale motions
4160	!================================================================
4161
4162	d_t_ec(:,:)=0.
4163	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4164	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
4165	& u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4166	& zmasse,exner,d_t_ec)
4167	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4168
4169	fname = 'after ener_conserv'
4170	varnamechk = 'paprs'
4171	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
4172	varnamechk = 'pplay'
4173	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
4174	varnamechk = 'pphi'
4175	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
4176	varnamechk = 't_seri'
4177	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
4178	varnamechk = 'u_seri'
4179	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
4180	varnamechk = 'v_seri'
4181	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
4182	varnamechk = 'q_seri'
4183	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
4184	varnamechk = 'd_t_ec'
4185	CALL check_var3D(fname, varnamechk, d_t_ec, klon, klev, largest, .FALSE.)
4186
4187	!IM
4188	IF (ip_ebil_phy.ge.1) THEN
4189	ztit='after physic'
4190	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
4191	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
4192	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
4193	!C Comme les tendances de la physique sont ajoute dans la dynamique,
4194	!C on devrait avoir que la variation d'entalpie par la dynamique
4195	!C est egale a la variation de la physique au pas de temps precedent.
4196	!C Donc la somme de ces 2 variations devrait etre nulle.
4197
4198	PRINT *,' LLuis sending: d_h_vcol: ',d_h_vcol,' d_qt ',d_qt,' d_ec ',d_ec
4199
4200	call diagphy(airephy,ztit,ip_ebil_phy &
4201	& , topsw, toplw, solsw, sollw, sens &
4202	& , evap, rain_fall, snow_fall, ztsol &
4203	& , d_h_vcol, d_qt, d_ec &
4204	& , fs_bound, fq_bound )
4205	!C
4206	d_h_vcol_phy=d_h_vcol
4207	!C
4208	END IF
4209	PRINT *,' Lluis 11 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
4210	fname = 'After everything, writting'
4211	varnamechk = 'paprs'
4212	CALL check_var3D(fname, varnamechk, paprs, klon, klev+1, largest, .FALSE.)
4213	varnamechk = 'pplay'
4214	CALL check_var3D(fname, varnamechk, pplay, klon, klev, largest, .FALSE.)
4215	varnamechk = 'pphi'
4216	CALL check_var3D(fname, varnamechk, pphi, klon, klev, largest, .FALSE.)
4217	varnamechk = 't_seri'
4218	CALL check_var3D(fname, varnamechk, t_seri, klon, klev, largest, .FALSE.)
4219	varnamechk = 'u_seri'
4220	CALL check_var3D(fname, varnamechk, u_seri, klon, klev, largest, .FALSE.)
4221	varnamechk = 'v_seri'
4222	CALL check_var3D(fname, varnamechk, v_seri, klon, klev, largest, .FALSE.)
4223	varnamechk = 'q_seri'
4224	CALL check_var3D(fname, varnamechk, q_seri, klon, klev, largest, .FALSE.)
4225	varnamechk = 'qv'
4226	CALL threscheck_var3D(fname, varnamechk, qx(:,:,1), klon, klev, 0.,-1, .FALSE.)
4227	varnamechk = 'q_seri'
4228	CALL threscheck_var3D(fname, varnamechk, q_seri, klon, klev, 0., -1, .FALSE.)
4229	varnamechk = 'pphi'
4230	CALL threscheck_var3D(fname, varnamechk, pphi, klon, klev, 0., -1, .FALSE.)
4231
4232	PRINT *,'Lluis Reaching the SORTIES point'
4233
4234	!C
4235	!c=======================================================================
4236	!c SORTIES
4237	!c=======================================================================
4238
4239	!IM Interpolation sur les niveaux de pression du NMC
4240	!c -------------------------------------------------
4241	!c
4242	#include "calcul_STDlev.h"
4243	!c
4244	!c slp sea level pressure
4245	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
4246	!c
4247	!ccc prw = eau precipitable
4248	DO i = 1, klon
4249	prw(i) = 0.
4250	DO k = 1, klev
4251	prw(i) = prw(i) + &
4252	& q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
4253	ENDDO
4254	ENDDO
4255	!c
4256	!IM initialisation + calculs divers diag AMIP2
4257	!c
4258	#include "calcul_divers.h"
4259	!c
4260	IF (type_trac == 'inca') THEN
4261	#ifdef INCA
4262	CALL VTe(VTphysiq)
4263	CALL VTb(VTinca)
4264
4265	CALL chemhook_end ( &
4266	& dtime, &
4267	& pplay, &
4268	& t_seri, &
4269	& tr_seri, &
4270	& nbtr, &
4271	& paprs, &
4272	& q_seri, &
4273	& airephy, &
4274	& pphi, &
4275	& pphis, &
4276	& zx_rh)
4277
4278	CALL VTe(VTinca)
4279	CALL VTb(VTphysiq)
4280	#endif
4281	END IF
4282
4283
4284	!c
4285	!c Convertir les incrementations en tendances
4286	!c
4287	IF (prt_level .GE.10) THEN
4288	print *,'Convertir les incrementations en tendances '
4289	ENDIF
4290	!c
4291	if (mydebug) then
4292	call writefield_phy('u_seri',u_seri,llm)
4293	call writefield_phy('v_seri',v_seri,llm)
4294	call writefield_phy('t_seri',t_seri,llm)
4295	call writefield_phy('q_seri',q_seri,llm)
4296	endif
4297
4298	DO k = 1, klev
4299	DO i = 1, klon
4300	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4301	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4302	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4303	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4304	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4305	ENDDO
4306	ENDDO
4307	!c
4308	IF (nqtot.GE.3) THEN
4309	DO iq = 3, nqtot
4310	DO k = 1, klev
4311	DO i = 1, klon
4312	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4313	ENDDO
4314	ENDDO
4315	ENDDO
4316	ENDIF
4317	!c
4318	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4319	!IM global posePB#include "write_bilKP_ins.h"
4320	!IM global posePB#include "write_bilKP_ave.h"
4321	!c
4322
4323	!c Sauvegarder les valeurs de t et q a la fin de la physique:
4324	!c
4325	DO k = 1, klev
4326	DO i = 1, klon
4327	u_ancien(i,k) = u_seri(i,k)
4328	v_ancien(i,k) = v_seri(i,k)
4329	t_ancien(i,k) = t_seri(i,k)
4330	q_ancien(i,k) = q_seri(i,k)
4331	ENDDO
4332	ENDDO
4333
4334	!!! RomP >>>
4335	IF (nqtot.GE.3) THEN
4336	DO iq = 3, nqtot
4337	DO k = 1, klev
4338	DO i = 1, klon
4339	tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
4340	ENDDO
4341	ENDDO
4342	ENDDO
4343	ENDIF
4344	!!! RomP <<<
4345	!==========================================================================
4346	! Sorties des tendances pour un point particulier
4347	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4348	! pour le debug
4349	! La valeur de igout est attribuee plus haut dans le programme
4350	!==========================================================================
4351
4352	if (prt_level.ge.1) then
4353	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4354	write(lunout,*) &
4355	& 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4356	write(lunout,*) &
4357	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4358	& pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4359	& pctsrf(igout,is_sic)
4360	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4361	do k=1,klev
4362	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4363	& d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4364	& d_t_eva(igout,k)
4365	enddo
4366	write(lunout,*) 'cool,heat'
4367	do k=1,klev
4368	write(lunout,*) cool(igout,k),heat(igout,k)
4369	enddo
4370
4371	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4372	do k=1,klev
4373	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4374	& d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4375	enddo
4376
4377	write(lunout,*) 'd_ps ',d_ps(igout)
4378	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4379	do k=1,klev,klon
4380	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4381	& d_qx(igout,k,1),d_qx(igout,k,2)
4382	enddo
4383	endif
4384
4385	!==========================================================================
4386
4387	!c============================================================
4388	!c Calcul de la temperature potentielle
4389	!c============================================================
4390	DO k = 1, klev
4391	DO i = 1, klon
4392	!cJYG/IM theta en debut du pas de temps
4393	!cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4394	!cJYG/IM theta en fin de pas de temps de physique
4395	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4396	!c thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
4397	!c fth_fonctions.F90 et parkind1.F90
4398	!c sinon thetal=theta
4399	!c thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4400	!c : ql_seri(i,k))
4401	thetal(i,k)=theta(i,k)
4402	ENDDO
4403	ENDDO
4404	!c
4405
4406	!c 22.03.04 BEG
4407	!c=============================================================
4408	!c Ecriture des sorties
4409	!c=============================================================
4410	#ifdef CPP_IOIPSL
4411
4412	!c Recupere des varibles calcule dans differents modules
4413	!c pour ecriture dans histxxx.nc
4414
4415	! Get some variables from module fonte_neige_mod
4416	!L. Fita, LMD. November 2013
4417	!! It is not working after removing starphy.nc and limit.nc?
4418	CALL fonte_neige_get_vars(pctsrf, &
4419	& zxfqcalving, zxfqfonte, zxffonte)
4420
4421
4422
4423	!c=============================================================
4424	! Separation entre thermiques et non thermiques dans les sorties
4425	! de fisrtilp
4426	!c=============================================================
4427
4428	if (iflag_thermals>=1) then
4429	d_t_lscth=0.
4430	d_t_lscst=0.
4431	d_q_lscth=0.
4432	d_q_lscst=0.
4433	do k=1,klev
4434	do i=1,klon
4435	if (ptconvth(i,k)) then
4436	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4437	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4438	else
4439	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4440	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4441	endif
4442	enddo
4443	enddo
4444
4445	do i=1,klon
4446	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4447	plul_th(i)=prfl(i,1)+psfl(i,1)
4448	enddo
4449	endif
4450
4451	PRINT *,' Lluis WRITING outputs! itap: ',itap,' itau_phy: ',itau_phy, &
4452	' itau_w: ',itau_w
4453	PRINT *,' Lluis writting outputs qsol: ',qsol(llp), &
4454	' ftsol: ',ftsol(llp,:)
4455
4456	#include "phys_output_write_new.h"
4457
4458
4459	#ifdef histISCCP
4460	#include "write_histISCCP.h"
4461	#endif
4462
4463	PRINT *,' Lluis WRITING histfiles!'
4464
4465	#ifdef histNMC
4466	#include "write_histhfNMC.h"
4467	#include "write_histdayNMC.h"
4468	#include "write_histmthNMC.h"
4469	#endif
4470
4471	#include "write_histday_seri.h"
4472
4473	#include "write_paramLMDZ_phy.h"
4474
4475	#endif
4476
4477	!c 22.03.04 END
4478	!c
4479	!c====================================================================
4480	!c Si c'est la fin, il faut conserver l'etat de redemarrage
4481	!c====================================================================
4482	!c
4483
4484	!c -----------------------------------------------------------------
4485	!c WSTATS: Saving statistics
4486	!c -----------------------------------------------------------------
4487	!c ("stats" stores and accumulates 8 key variables in file "stats.nc"
4488	!c which can later be used to make the statistic files of the run:
4489	!c "stats") only possible in 3D runs !
4490
4491
4492	IF (callstats) THEN
4493
4494	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
4495	& ,2,paprs(:,1))
4496	call wstats(klon,o_tsol%name,"Surface temperature","K", &
4497	& 2,zxtsol)
4498	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4499	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
4500	& "kg/(s*m2)",2,zx_tmp_fi2d)
4501	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4502	call wstats(klon,o_plul%name,"Large-scale Precip", &
4503	& "kg/(s*m2)",2,zx_tmp_fi2d)
4504	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4505	call wstats(klon,o_pluc%name,"Convective Precip", &
4506	& "kg/(s*m2)",2,zx_tmp_fi2d)
4507	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
4508	& "W/m2",2,solsw)
4509	call wstats(klon,o_soll%name,"IR rad. at surf.", &
4510	& "W/m2",2,sollw)
4511	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4512	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
4513	& "W/m2",2,zx_tmp_fi2d)
4514
4515
4516
4517	call wstats(klon,o_temp%name,"Air temperature","K", &
4518	& 3,t_seri)
4519	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
4520	& 3,u_seri)
4521	call wstats(klon,o_vitv%name,"Meridional wind", &
4522	& "m.s-1",3,v_seri)
4523	call wstats(klon,o_vitw%name,"Vertical wind", &
4524	& "m.s-1",3,omega)
4525	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
4526	& 3,q_seri)
4527
4528
4529
4530	IF(lafin) THEN
4531	write (,) "Writing stats..."
4532	call mkstats(ierr)
4533	ENDIF
4534
4535	ENDIF !if callstats
4536
4537	IF (lafin) THEN
4538	itau_phy = itau_phy + itap
4539	CALL phyredem ("restartphy.nc")
4540	! open(97,form="unformatted",file="finbin")
4541	! write(97) u_seri,v_seri,t_seri,q_seri
4542	! close(97)
4543	!$OMP MASTER
4544	if (read_climoz >= 1) then
4545	if (is_mpi_root) then
4546	call nf95_close(ncid_climoz)
4547	end if
4548	deallocate(press_climoz) ! pointer
4549	end if
4550	!$OMP END MASTER
4551	ENDIF
4552
4553	! first=.false.
4554
4555	! Lluis
4556	PRINT *,' Lluis: ',klev,' UBOUNDS: ',UBOUND(t_seri), UBOUND(u_seri), &
4557	UBOUND(d_q_con), UBOUND(d_t_con)
4558	PRINT *,' Lluis llp ',llp,' itap: ',itap,' zlev t_seri u_seri d_q_con d_t_con_____'
4559	DO i=1,klev
4560	PRINT *,i,t_seri(llp,i), u_seri(llp,i), d_q_con(llp,i), d_t_con(llp,i)
4561	END DO
4562
4563
4564	RETURN
4565	END SUBROUTINE physiq
4566
4567	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4568	IMPLICIT none
4569	!c
4570	!c Calculer et imprimer l'eau totale. A utiliser pour verifier
4571	!c la conservation de l'eau
4572	!c
4573	#include "YOMCST.h"
4574	INTEGER klon,klev
4575	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4576	REAL aire(klon)
4577	REAL qtotal, zx, qcheck
4578	INTEGER i, k
4579	!c
4580	zx = 0.0
4581	DO i = 1, klon
4582	zx = zx + aire(i)
4583	ENDDO
4584	qtotal = 0.0
4585	DO k = 1, klev
4586	DO i = 1, klon
4587	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4588	& *(paprs(i,k)-paprs(i,k+1))/RG
4589	ENDDO
4590	ENDDO
4591	!c
4592	qcheck = qtotal/zx
4593	!c
4594	RETURN
4595	END
4596
4597	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4598	IMPLICIT none
4599	!c
4600	!c Tranformer une variable de la grille physique a
4601	!c la grille d'ecriture
4602	!c
4603	INTEGER nfield,nlon,iim,jjmp1, jjm
4604	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4605	!c
4606	INTEGER i, n, ig
4607	!c
4608	jjm = jjmp1 - 1
4609	DO n = 1, nfield
4610	DO i=1,iim
4611	ecrit(i,n) = fi(1,n)
4612	ecrit(i+jjm*iim,n) = fi(nlon,n)
4613	ENDDO
4614	DO ig = 1, nlon - 2
4615	ecrit(iim+ig,n) = fi(1+ig,n)
4616	ENDDO
4617	ENDDO
4618	RETURN
4619	END SUBROUTINE gr_fi_ecrit
4620
4621	SUBROUTINE check_var(funcn, varn, var, sizev, bigvalue, stoprun)
4622	! Subroutine to check the consistency of a variable
4623	! * NaN value: by definition is variable /= variable
4624	! * bigvalue: threshold for the variable
4625
4626	IMPLICIT NONE
4627
4628	#include "dimensions.h"
4629
4630	INTEGER, INTENT(IN) :: sizev
4631	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
4632	REAL, DIMENSION(sizev), INTENT(IN) :: var
4633	REAL, INTENT(IN) :: bigvalue
4634	LOGICAL, INTENT(IN) :: stoprun
4635
4636	! Local
4637	INTEGER :: i, wrongi, xpt, ypt
4638	CHARACTER(LEN=50) :: errmsg
4639	LOGICAL :: found
4640	REAL, DIMENSION(sizev) :: wrongvalues
4641	INTEGER, DIMENSION(sizev) :: wronggridpt
4642
4643	!!!!!!! Variables
4644	! funcn: at which functino of part of the program variable is checked
4645	! varn: name of the variable
4646	! var: variable to check
4647	! sizev: size of the variable
4648	! bigvalue: biggest attenaible value for the variable
4649	! stoprun: Should the run stop if it founds a problem?
4650
4651	errmsg = 'ERROR -- error -- ERROR -- error'
4652
4653	found = .FALSE.
4654	wrongi = 0
4655	DO i=1,sizev
4656	IF (var(i) /= var(i) .OR. ABS(var(i)) > bigvalue ) THEN
4657	IF (wrongi == 0) found = .TRUE.
4658	wrongi = wrongi + 1
4659	wrongvalues(wrongi) = var(i)
4660	wronggridpt(wrongi) = i
4661	END IF
4662	END DO
4663
4664	IF (found) THEN
4665	PRINT *,TRIM(errmsg)
4666	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
4667	"' is wrong in Nvalues= ",wrongi,' at i (x, y) value___'
4668	DO i=1,wrongi
4669	ypt = INT(wronggridpt(i)/wiim) + 1
4670	xpt = wronggridpt(i) - (ypt-1)*wiim
4671	PRINT *,wronggridpt(i), '(',xpt,', ',ypt,')', wrongvalues(i)
4672	END DO
4673	IF (stoprun) THEN
4674	STOP
4675	END IF
4676	END IF
4677
4678	RETURN
4679
4680	END SUBROUTINE check_var
4681
4682	SUBROUTINE check_var3D(funcn, varn, var, sizev, zsize, bigvalue, stoprun)
4683	! Subroutine to check the consistency of a 3D LMDSZ - variable (klon, klev) !
4684	! * NaN value: by definition is variable /= variable
4685	! * bigvalue: threshold for the variable
4686
4687	IMPLICIT NONE
4688
4689	#include "dimensions.h"
4690
4691	INTEGER, INTENT(IN) :: sizev, zsize
4692	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
4693	REAL, DIMENSION(sizev,zsize), INTENT(IN) :: var
4694	REAL, INTENT(IN) :: bigvalue
4695	LOGICAL, INTENT(IN) :: stoprun
4696
4697	! Local
4698	INTEGER :: i, k, wrongi, xpt, ypt
4699	CHARACTER(LEN=50) :: errmsg
4700	LOGICAL :: found
4701	REAL, DIMENSION(sizev*zsize) :: wrongvalues
4702	INTEGER, DIMENSION(sizev*zsize,2) :: wronggridpt
4703
4704	!!!!!!! Variables
4705	! funcn: at which functino of part of the program variable is checked
4706	! varn: name of the variable
4707	! var: variable to check
4708	! sizev: size of the variable
4709	! zsize: vertical size of the variable
4710	! bigvalue: biggest attenaible value for the variable
4711	! stoprun: Should the run stop if it founds a problem?
4712
4713	errmsg = 'ERROR -- error -- ERROR -- error'
4714
4715	found = .FALSE.
4716	wrongi = 0
4717	DO i=1,sizev
4718	DO k=1,zsize
4719	IF (var(i,k) /= var(i,k) .OR. ABS(var(i,k)) > bigvalue ) THEN
4720	IF (wrongi == 0) found = .TRUE.
4721	wrongi = wrongi + 1
4722	wrongvalues(wrongi) = var(i,k)
4723	wronggridpt(wrongi,1) = i
4724	wronggridpt(wrongi,2) = k
4725	END IF
4726	END DO
4727	END DO
4728
4729	IF (found) THEN
4730	PRINT *,TRIM(errmsg)
4731	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
4732	"' is wrong in Nvalues= ",wrongi,' at i (x,y) k value___'
4733	DO i=1,wrongi
4734	ypt = INT(wronggridpt(i,1)/wiim) + 1
4735	xpt = wronggridpt(i,1) - (ypt-1)*wiim
4736	PRINT *,wronggridpt(i,1), '(',xpt,', ',ypt,')', wronggridpt(i,2), wrongvalues(i)
4737	END DO
4738	IF (stoprun) THEN
4739	STOP
4740	END IF
4741	END IF
4742
4743	RETURN
4744
4745	END SUBROUTINE check_var3D
4746
4747	SUBROUTINE threscheck_var3D(funcn, varn, var, sizev, zsize, thresvalue, kthres, &
4748	stoprun)
4749	! Subroutine to check the consistency of a 3D LMDSZ - variable (klon, klev) !
4750	! * thresvalue: threshold for the variable
4751
4752	IMPLICIT NONE
4753
4754	#include "dimensions.h"
4755
4756	INTEGER, INTENT(IN) :: sizev, zsize, kthres
4757	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
4758	REAL, DIMENSION(sizev,zsize), INTENT(IN) :: var
4759	REAL, INTENT(IN) :: thresvalue
4760	LOGICAL, INTENT(IN) :: stoprun
4761
4762	! Local
4763	INTEGER :: i, k, wrongi, xpt, ypt
4764	CHARACTER(LEN=50) :: errmsg
4765	LOGICAL :: found
4766	REAL, DIMENSION(sizev*zsize) :: wrongvalues
4767	INTEGER, DIMENSION(sizev*zsize,2) :: wronggridpt
4768
4769	!!!!!!! Variables
4770	! funcn: at which functino of part of the program variable is checked
4771	! varn: name of the variable
4772	! var: variable to check
4773	! sizev: size of the variable
4774	! zsize: vertical size of the variable
4775	! thresvalue: threshold attenaible value for the variable
4776	! kthres: kind of threshold= -1, below wrong, 1: above wrong
4777	! stoprun: Should the run stop if it founds a problem?
4778
4779	errmsg = 'ERROR -- error -- ERROR -- error'
4780
4781	found = .FALSE.
4782	wrongi = 0
4783	DO i=1,sizev
4784	DO k=1,zsize
4785	IF (kthres == -1 .AND. var(i,k) < thresvalue ) THEN
4786	IF (wrongi == 0) found = .TRUE.
4787	wrongi = wrongi + 1
4788	wrongvalues(wrongi) = var(i,k)
4789	wronggridpt(wrongi,1) = i
4790	wronggridpt(wrongi,2) = k
4791	ELSE IF (kthres == 1 .AND. var(i,k) > thresvalue ) THEN
4792	IF (wrongi == 0) found = .TRUE.
4793	wrongi = wrongi + 1
4794	wrongvalues(wrongi) = var(i,k)
4795	wronggridpt(wrongi,1) = i
4796	wronggridpt(wrongi,2) = k
4797	END IF
4798	END DO
4799	END DO
4800
4801	IF (found) THEN
4802	PRINT *,TRIM(errmsg)
4803	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
4804	"' is wrong (threshold: ',thresvalue,') in Nvalues= ",wrongi, &
4805	' at i (x,y) k value___'
4806	DO i=1,wrongi
4807	ypt = INT(wronggridpt(i,1)/wiim) + 1
4808	xpt = wronggridpt(i,1) - (ypt-1)*wiim
4809	PRINT *,wronggridpt(i,1), '(',xpt,', ',ypt,')', wronggridpt(i,2), wrongvalues(i)
4810	END DO
4811	IF (stoprun) THEN
4812	STOP
4813	END IF
4814	END IF
4815
4816	RETURN
4817
4818	END SUBROUTINE threscheck_var3D
4819

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

Download in other formats:

Original Format