Context Navigation

physiq.F90 @ 2262

Last change on this file since 2262 was 2262, checked in by jyg, 10 years ago

Correction of a bug concerning the number 'nqo' of
water phases transported by the dynamic : the
default value (= 2, corresponding to vapour and
liquid phases) was still explicitely present in
various places.

Modified files:

infotrac.F90,
physiq.F90

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: 138.0 KB

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

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

Download in other formats:

Original Format