Context Navigation

physiq_mod.F90 @ 2569

Last change on this file since 2569 was 2569, checked in by Laurent Fairhead, 8 years ago
Pour retrouver 1+1=2 avec iflag_albedo=1 LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 149.3 KB

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

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

Download in other formats:

Original Format