Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq.F90

Last change on this file was 164, checked in by lfita, 10 years ago
Checking more properly 'q_seri'
File size: 187.4 KB

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

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

Download in other formats:

Original Format