Context Navigation

physiq.F90 @ 113

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

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

Download in other formats:

Original Format