Context Navigation

physiq_mod.F @ 3726

Last change on this file since 3726 was 3726, checked in by emillour, 3 months ago
Mars PCM: Turn "callkeys.h" into module "callkeys_mod.F90" EM
File size: 163.1 KB

Rev	Line
[1549]	1	MODULE physiq_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[234]	7	SUBROUTINE physiq(
[226]	8	$ ngrid,nlayer,nq
	9	$ ,firstcall,lastcall
	10	$ ,pday,ptime,ptimestep
	11	$ ,pplev,pplay,pphi
	12	$ ,pu,pv,pt,pq
[1312]	13	$ ,flxw
[1576]	14	$ ,pdu,pdv,pdt,pdq,pdpsrf)
[42]	15
[2162]	16	use watercloud_mod, only: watercloud, zdqcloud, zdqscloud
[2164]	17	use calchim_mod, only: calchim, ichemistry, zdqchim, zdqschim
[1996]	18	use watersat_mod, only: watersat
[3207]	19	use co2condens_mod, only: co2condens, CO2cond_ps
[2562]	20	use co2cloud_mod, only: co2cloud
[1974]	21	use callradite_mod, only: callradite
	22	use callsedim_mod, only: callsedim
	23	use rocketduststorm_mod, only: rocketduststorm, dustliftday
	24	use calcstormfract_mod, only: calcstormfract
[2628]	25	use topmons_mod, only: topmons,topmons_setup
[3012]	26	use nltecool_mod, only: nltecool
	27	use nlte_tcool_mod, only: nlte_tcool
[3016]	28	use blendrad_mod, only: blendrad
	29	use nlthermeq_mod, only: nlthermeq
[3015]	30	use thermosphere_mod, only: thermosphere
[3466]	31	use param_read_e107_mod, only: param_read_e107
[1617]	32	use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
[1036]	33	& igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
[2312]	34	& igcm_hdo_vap, igcm_hdo_ice,
[1036]	35	& igcm_ccn_mass, igcm_ccn_number,
[1617]	36	& igcm_ccnco2_mass, igcm_ccnco2_number,
[2562]	37	& igcm_ccnco2_h2o_mass_ice,
	38	& igcm_ccnco2_h2o_mass_ccn,
	39	& igcm_ccnco2_h2o_number,
[2589]	40	& igcm_ccnco2_meteor_mass,
	41	& igcm_ccnco2_meteor_number,
[1036]	42	& igcm_dust_mass, igcm_dust_number, igcm_h2o2,
[1660]	43	& nuice_ref, rho_ice, rho_dust, ref_r0,
[1974]	44	& igcm_he, igcm_stormdust_mass,
[2199]	45	& igcm_stormdust_number, igcm_topdust_mass,
[2324]	46	& igcm_topdust_number,
	47	& qperemin
[2942]	48	use comsoil_h, only: inertiedat, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start
[2285]	49	& tsoil, nsoilmx,!number of subsurface layers
[3113]	50	& mlayer,layer, ! soil mid layer depths
	51	& nqsoil,qsoil ! adsorption
[1974]	52	use geometry_mod, only: longitude, latitude, cell_area,
[3078]	53	& cell_area_for_lonlat_outputs,longitude_deg
[1541]	54	use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
[3203]	55	use surfdat_h, only: phisfi, albedodat, z0, albedo_h2o_cap,
	56	& albedo_h2o_frost, frost_albedo_threshold,
	57	& frost_metam_threshold, tsurf, emis, capcal,
	58	& fluxgrd, qsurf, watercap, watercaptag,
[3130]	59	& perennial_co2ice
[3203]	60	use comsaison_h, only: dist_sol, declin, zls,
[2597]	61	& mu0, fract, local_time
[3094]	62	use solarlong_mod, only: solarlong
[3006]	63	use nirdata_mod, only: NIR_leedat
	64	use nirco2abs_mod, only: nirco2abs
[3466]	65	use surfacearea_mod, only: surfacearea
[3183]	66	use slope_mod, only: theta_sl, psi_sl, getslopes, param_slope
[3203]	67	use conc_mod, only: init_r_cp_mu, update_r_cp_mu_ak, rnew,
[3185]	68	& cpnew, mmean
[3369]	69	use time_phylmdz_mod, only: steps_per_sol
	70	use time_phylmdz_mod, only: iphysiq, ecritstart, daysec
[2409]	71	use dimradmars_mod, only: aerosol, totcloudfrac,
[1246]	72	& dtrad, fluxrad_sky, fluxrad, albedo,
[3203]	73	& naerkind, iaer_dust_doubleq,
[2678]	74	& iaer_stormdust_doubleq, iaer_h2o_ice,
	75	& flux_1AU
[2409]	76	use dust_param_mod, only: doubleq, lifting, callddevil,
[2417]	77	& tauscaling, odpref, dustbin,
[2442]	78	& dustscaling_mode, dust_rad_adjust,
[2643]	79	& freedust, reff_driven_IRtoVIS_scenario
[3292]	80	use dustdevil_mod, only: dustdevil
[3203]	81	use turb_mod, only: q2, wstar, ustar, sensibFlux,
[1236]	82	& zmax_th, hfmax_th, turb_resolved
[1229]	83	use planete_h, only: aphelie, periheli, year_day, peri_day,
	84	& obliquit
[3095]	85	use planete_h, only: iniorbit
[3203]	86	USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
[1912]	87	USE calldrag_noro_mod, ONLY: calldrag_noro
[3262]	88	USE vdifc_mod, ONLY: vdifc
[3203]	89	use param_v4_h, only: nreact,n_avog,
[1266]	90	& fill_data_thermos, allocate_param_thermos
	91	use iono_h, only: allocate_param_iono
[1974]	92	use compute_dtau_mod, only: compute_dtau
[2149]	93	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
[3144]	94	use nonoro_gwd_mix_mod, only: nonoro_gwd_mix, calljliu_gwimix
[2551]	95	use check_fields_mod, only: check_physics_fields
[3142]	96	use surfini_mod, only: surfini
[1112]	97	#ifdef MESOSCALE
	98	use comsoil_h, only: mlayer,layer
	99	use surfdat_h, only: z0_default
[1236]	100	use comm_wrf
[1212]	101	#else
[2223]	102	USE planetwide_mod, ONLY: planetwide_maxval, planetwide_minval,
	103	& planetwide_sumval
[1212]	104	use phyredem, only: physdem0, physdem1
[2545]	105	use phyetat0_mod, only: phyetat0, tab_cntrl_mod
[2559]	106	use wstats_mod, only: callstats, wstats, mkstats
[1212]	107	use eofdump_mod, only: eofdump
[2223]	108	USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
	109	USE mod_phys_lmdz_omp_data, ONLY: is_omp_master
[2511]	110	USE time_phylmdz_mod, ONLY: day_end
[1112]	111	#endif
	112
[3203]	113	#ifdef CPP_XIOS
[2223]	114	use xios_output_mod, only: initialize_xios_output,
	115	& update_xios_timestep,
	116	& send_xios_field
	117	use wxios, only: wxios_context_init, xios_context_finalize
	118	#endif
[3203]	119	USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured,
[3078]	120	& regular_lonlat
[2551]	121	use ioipsl_getin_p_mod, only: getin_p
[2896]	122	use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
	123	& subslope_dist,iflat,sky_slope,
	124	& major_slope,compute_meshgridavg,
	125	& ini_comslope_h
[2932]	126	use write_output_mod, only: write_output
[3151]	127	use pbl_parameters_mod, only: pbl_parameters
[3167]	128	use lmdz_atke_turbulence_ini, only : atke_ini
[3230]	129	use waterice_tifeedback_mod, only : waterice_tifeedback
[3726]	130	use callkeys_mod, only: calladj, calltherm, callatke, calldifv
	131	use callkeys_mod, only: callrichsl, tke_heat_flux
	132	use callkeys_mod, only: calllott, calllott_nonoro, calleofdump
	133	use callkeys_mod, only: callrad, callnlte, callnirco2, nircorr
	134	use callkeys_mod, only: diurnal, season, iradia, nltemodel
	135	use callkeys_mod, only: water, activice, microphys, CLFvarying
	136	use callkeys_mod, only: hdo, co2clouds, co2useh2o, meteo_flux
	137	use callkeys_mod, only: callsoil, callslope, callcond
	138	use callkeys_mod, only: tituscap, surfaceice_tifeedback
	139	use callkeys_mod, only: refill_watercap, poreice_tifeedback
	140	use callkeys_mod, only: cst_cap_albedo
	141	use callkeys_mod, only: rdstorm, dustinjection
	142	use callkeys_mod, only: topflows, dustiropacity
	143	use callkeys_mod, only: sedimentation, scavenging
	144	use callkeys_mod, only: photochem, callthermos
	145	use callkeys_mod, only: startphy_file
	146
[42]	147	IMPLICIT NONE
	148	c=======================================================================
	149	c
	150	c subject:
	151	c --------
	152	c
[3203]	153	c Organisation of the physical parametrisations of the LMD
[42]	154	c martian atmospheric general circulation model.
	155	c
	156	c The GCM can be run without or with tracer transport
	157	c depending on the value of Logical "tracer" in file "callphys.def"
	158	c Tracers may be water vapor, ice OR chemical species OR dust particles
	159	c
	160	c SEE comments in initracer.F about numbering of tracer species...
	161	c
	162	c It includes:
	163	c
	164	c 1. Initialization:
	165	c 1.1 First call initializations
	166	c 1.2 Initialization for every call to physiq
	167	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	168	c 2. Compute radiative transfer tendencies
	169	c (longwave and shortwave) for CO2 and aerosols.
	170	c 3. Gravity wave and subgrid scale topography drag :
	171	c 4. Vertical diffusion (turbulent mixing):
	172	c 5. Convective adjustment
	173	c 6. Condensation and sublimation of carbon dioxide.
	174	c 7. TRACERS :
[1617]	175	c 7a. water, water ice, co2 ice (clouds)
[42]	176	c 7b. call for photochemistry when tracers are chemical species
	177	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
	178	c 7d. updates (CO2 pressure variations, surface budget)
	179	c 8. Contribution to tendencies due to thermosphere
	180	c 9. Surface and sub-surface temperature calculations
	181	c 10. Write outputs :
	182	c - "startfi", "histfi" (if it's time)
	183	c - Saving statistics (if "callstats = .true.")
	184	c - Dumping eof (if "calleofdump = .true.")
[3203]	185	c - Output any needed variables in "diagfi"
[42]	186	c 11. Diagnostic: mass conservation of tracers
[3203]	187	c
	188	c author:
	189	c -------
	190	c Frederic Hourdin 15/10/93
	191	c Francois Forget 1994
	192	c Christophe Hourdin 02/1997
[42]	193	c Subroutine completly rewritten by F.Forget (01/2000)
	194	c Introduction of the photochemical module: S. Lebonnois (11/2002)
	195	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
	196	c Water ice clouds: Franck Montmessin (update 06/2003)
	197	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	198	c Nb: See callradite.F for more information.
[234]	199	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
[414]	200	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
[1617]	201	c
	202	c 10/16 J. Audouard: modifications for CO2 clouds scheme
	203
[42]	204	c arguments:
	205	c ----------
	206	c
	207	c input:
	208	c ------
	209	c ecri period (in dynamical timestep) to write output
	210	c ngrid Size of the horizontal grid.
	211	c All internal loops are performed on that grid.
	212	c nlayer Number of vertical layers.
	213	c nq Number of advected fields
	214	c firstcall True at the first call
	215	c lastcall True at the last call
	216	c pday Number of days counted from the North. Spring
	217	c equinoxe.
	218	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	219	c ptimestep timestep (s)
	220	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	221	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	222	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	223	c pu(ngrid,nlayer) u component of the wind (ms-1)
	224	c pv(ngrid,nlayer) v component of the wind (ms-1)
	225	c pt(ngrid,nlayer) Temperature (K)
	226	c pq(ngrid,nlayer,nq) Advected fields
[3203]	227	c pudyn(ngrid,nlayer) \|
[330]	228	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
	229	c ptdyn(ngrid,nlayer) \| corresponding variables
	230	c pqdyn(ngrid,nlayer,nq) \|
[1312]	231	c flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
[42]	232	c
	233	c output:
	234	c -------
	235	c
[1047]	236	c pdu(ngrid,nlayer) \|
	237	c pdv(ngrid,nlayer) \| Temporal derivative of the corresponding
	238	c pdt(ngrid,nlayer) \| variables due to physical processes.
	239	c pdq(ngrid,nlayer,nq) \|
	240	c pdpsrf(ngrid) \|
[42]	241
	242	c
	243	c=======================================================================
	244	c
	245	c 0. Declarations :
	246	c ------------------
	247
[1922]	248	include "netcdf.inc"
[42]	249
	250	c Arguments :
	251	c -----------
	252
	253	c inputs:
	254	c -------
[1312]	255	INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
	256	INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
	257	INTEGER,INTENT(in) :: nq ! number of tracers
	258	LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
	259	LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
	260	REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
	261	REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
	262	REAL,INTENT(in) :: ptimestep ! physics timestep (s)
	263	REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	264	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	265	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
	266	REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
	267	REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
	268	REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
	269	REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
[3203]	270	REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s) at lower boundary of layer
[226]	271
[42]	272	c outputs:
	273	c --------
	274	c physical tendencies
[1312]	275	REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
	276	REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
	277	REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
	278	REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
	279	REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
[42]	280
	281	c Local saved variables:
	282	c ----------------------
[1233]	283	INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
[3042]	284	INTEGER,SAVE :: icount ! Counter of calls to physiq during the run.
[2578]	285	REAL,SAVE :: time_phys
[1974]	286
[2578]	287	!$OMP THREADPRIVATE(day_ini,icount,time_phys)
	288
[1375]	289	#ifdef DUSTSTORM
	290	REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
	291	#endif
[3203]	292
[42]	293	c Variables used by the water ice microphysical scheme:
[1047]	294	REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m)
	295	REAL nuice(ngrid,nlayer) ! Estimated effective variance
[42]	296	! of the size distribution
[1047]	297	real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
	298	real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3)
[2447]	299	real rsedcloudco2(ngrid,nlayer) ! CO2 Cloud sedimentation radius
	300	real rhocloudco2(ngrid,nlayer) ! CO2 Cloud density (kg.m-3)
[3203]	301	real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the
[2447]	302	! size distribution
[2942]	303	REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface
[1047]	304	! thermal inertia (J.s-1/2.m-2.K-1)
[3230]	305	! (used only when tifeedback surface or pore =.true.)
[1617]	306	c Variables used by the CO2 clouds microphysical scheme:
[1934]	307	DOUBLE PRECISION riceco2(ngrid,nlayer) ! co2 ice geometric mean radius (m)
[1617]	308	real zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
[2562]	309	real zdqssed_ccn(ngrid,nq) ! CCN flux at the surface (kg.m-2.s-1)
[3203]	310	real, dimension(ngrid,nlayer) :: zcondicea_co2microp
[1495]	311	c Variables used by the photochemistry
	312	REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
	313	REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry)
[234]	314	c Variables used by the slope model
[3203]	315	REAL sl_lct, sl_lat
[234]	316	REAL sl_tau, sl_alb, sl_the, sl_psi
	317	REAL sl_fl0, sl_flu
	318	REAL sl_ra, sl_di0
	319	REAL sky
[2685]	320	REAL fluxsurf_dir_dn_sw(ngrid) ! Incident direct solar flux on Mars at surface (W.m-2)
[234]	321
[1047]	322	REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
[42]	323
	324	c Local variables :
	325	c -----------------
	326	REAL CBRT
	327	EXTERNAL CBRT
	328
[3203]	329	! CHARACTER*80 fichier
	330	INTEGER l,ig,ierr,igout,iq,isoil
[42]	331
[2900]	332	REAL fluxsurf_lw(ngrid,nslope) !incident LW (IR) surface flux (W.m-2)
	333	REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2)
[2685]	334	REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2)
[1047]	335	REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2)
[2685]	336	REAL fluxtop_dn_sw(ngrid,2) ! Incoming SW (solar) flux from space (W.m-2)
	337	REAL fluxtop_up_sw(ngrid,2) ! Outgoing SW (solar) flux to space (W.m-2)
[2415]	338	REAL tau_pref_scenario(ngrid) ! prescribed dust column visible opacity
	339	! at odpref
[2643]	340	REAL IRtoVIScoef(ngrid) ! conversion coefficient to apply on
	341	! scenario absorption IR (9.3um) CDOD
	342	! = tau_pref_gcm_VIS / tau_pref_gcm_IR
[2415]	343	REAL tau_pref_gcm(ngrid) ! dust column visible opacity at odpref in the GCM
[1974]	344	c rocket dust storm
	345	REAL totstormfract(ngrid) ! fraction of the mesh where the dust storm is contained
[3203]	346	logical clearatm ! clearatm used to calculate twice the radiative
	347	! transfer when rdstorm is active :
[1974]	348	! - in a mesh with stormdust and background dust (false)
	349	! - in a mesh with background dust only (true)
[2628]	350	c entrainment by mountain top dust flows
[3203]	351	logical nohmons ! nohmons used to calculate twice the radiative
	352	! transfer when topflows is active :
[2199]	353	! - in a mesh with topdust and background dust (false)
	354	! - in a mesh with background dust only (true)
[3203]	355
[1047]	356	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
[1246]	357	! AS: TBD: this one should be in a module !
[42]	358	REAL zday ! date (time since Ls=0, in martian days)
[1047]	359	REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers
	360	REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
[3157]	361	REAL gz(ngrid,nlayer) ! variation of g with altitude from aeroid surface
[1036]	362	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
[42]	363
	364	c Tendancies due to various processes:
[2900]	365	REAL dqsurf(ngrid,nq,nslope) ! tendency for tracers on surface (Kg/m2/s)
[1047]	366	REAL zdtlw(ngrid,nlayer) ! (K/s)
	367	REAL zdtsw(ngrid,nlayer) ! (K/s)
[1974]	368	REAL pdqrds(ngrid,nlayer,nq) ! tendency for dust after rocketduststorm
	369
[1047]	370	REAL zdtnirco2(ngrid,nlayer) ! (K/s)
	371	REAL zdtnlte(ngrid,nlayer) ! (K/s)
[2900]	372	REAL zdtsurf(ngrid,nslope) ! (K/s)
[1617]	373	REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
[1047]	374	REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
[2900]	375	REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope) ! (K/s)
[1047]	376	REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
	377	REAL zdhadj(ngrid,nlayer) ! (K/s)
	378	REAL zdtgw(ngrid,nlayer) ! (K/s)
	379	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
[2900]	380	REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope)
[1047]	381	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
[42]	382
[2900]	383	REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope)
[1047]	384	REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
	385	REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
	386	REAL zdqadj(ngrid,nlayer,nq)
	387	REAL zdqc(ngrid,nlayer,nq)
[2162]	388	REAL zdqcloudco2(ngrid,nlayer,nq)
[2900]	389	REAL zdqsc(ngrid,nq,nslope)
[42]	390
[1047]	391	REAL zdteuv(ngrid,nlayer) ! (K/s)
	392	REAL zdtconduc(ngrid,nlayer) ! (K/s)
	393	REAL zdumolvis(ngrid,nlayer)
	394	REAL zdvmolvis(ngrid,nlayer)
	395	real zdqmoldiff(ngrid,nlayer,nq)
[2467]	396	real*8 PhiEscH,PhiEscH2,PhiEscD
[42]	397
[2900]	398	REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope) ! (kg/m-2)
[2260]	399
[42]	400	c Local variable for local intermediate calcul:
[2900]	401	REAL zflubid(ngrid,nslope)
[1047]	402	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
	403	REAL zdum1(ngrid,nlayer)
	404	REAL zdum2(ngrid,nlayer)
[42]	405	REAL ztim1,ztim2,ztim3, z1,z2
	406	REAL ztime_fin
[1047]	407	REAL zdh(ngrid,nlayer)
[1313]	408	REAL zh(ngrid,nlayer) ! potential temperature (K)
[1312]	409	REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
[42]	410	INTEGER length
	411	PARAMETER (length=100)
[3157]	412	REAL tlaymean ! temporary value of mean layer temperature for zzlay
[42]	413
[1974]	414	c Variables for the total dust for diagnostics
	415	REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust
	416
[42]	417	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	418	c -----------------------------------------------------------------------------
[1047]	419	REAL ps(ngrid), zt(ngrid,nlayer)
	420	REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
	421	REAL zq(ngrid,nlayer,nq)
[1974]	422
[2685]	423	REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid)
[2900]	424	REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid)
[42]	425	character*2 str2
[1036]	426	! character*5 str5
[1047]	427	real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
	428	real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
[1974]	429	real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m)
[2199]	430	real rtopdust(ngrid,nlayer) ! topdust geometric mean radius (m)
[42]	431	integer igmin, lmin
	432
[883]	433	! pplev and pplay are dynamical inputs and must not be modified in the physics.
	434	! instead, use zplay and zplev :
[3203]	435	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
[1036]	436	! REAL zstress(ngrid),cd
[1047]	437	real rho(ngrid,nlayer) ! density
	438	real vmr(ngrid,nlayer) ! volume mixing ratio
	439	real rhopart(ngrid,nlayer) ! number density of a given species
	440	real colden(ngrid,nq) ! vertical column of tracers
[1464]	441	real mass(nq) ! global mass of tracers (g)
[1047]	442	REAL mtot(ngrid) ! Total mass of water vapor (kg/m2)
[1974]	443	REAL mstormdtot(ngrid) ! Total mass of stormdust tracer (kg/m2)
[2362]	444	REAL mdusttot(ngrid) ! Total mass of dust tracer (kg/m2)
[1047]	445	REAL icetot(ngrid) ! Total mass of water ice (kg/m2)
[2551]	446	REAL mtotco2(ngrid) ! Total mass of co2, including ice at the surface (kg/m2)
	447	REAL vaptotco2(ngrid) ! Total mass of co2 vapor (kg/m2)
	448	REAL icetotco2(ngrid) ! Total mass of co2 ice (kg/m2)
[1047]	449	REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2)
	450	REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2)
	451	REAL rave(ngrid) ! Mean water ice effective radius (m)
	452	REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1
	453	REAL tauTES(ngrid) ! column optical depth at 825 cm-1
[42]	454	REAL Qabsice ! Water ice absorption coefficient
[1047]	455	REAL taucloudtes(ngrid) ! Cloud opacity at infrared
[520]	456	! reference wavelength using
	457	! Qabs instead of Qext
	458	! (direct comparison with TES)
[2312]	459	REAL mtotD(ngrid) ! Total mass of HDO vapor (kg/m2)
	460	REAL icetotD(ngrid) ! Total mass of HDO ice (kg/m2)
	461	REAL DoH_vap(ngrid,nlayer) !D/H ratio
	462	REAL DoH_ice(ngrid,nlayer) !D/H ratio
	463	REAL DoH_surf(ngrid) !D/H ratio surface
[2362]	464
[1047]	465	REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
	466	REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
	467	REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
	468	REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
	469	REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg)
	470	REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg)
[1974]	471	c definition tendancies of stormdust tracers
	472	REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s)
	473	REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s)
	474	REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg)
	475	REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg)
[2414]	476	REAL wspeed(ngrid,nlayer+1) ! vertical velocity stormdust tracer
	477	REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer
[2413]	478	REAL dsodust(ngrid,nlayer) ! density scaled opacity for background dust
	479	REAL dsords(ngrid,nlayer) ! density scaled opacity for stormdust
	480	REAL dsotop(ngrid,nlayer) ! density scaled opacity for topdust
[1974]	481
[411]	482	c Test 1d/3d scavenging
[1047]	483	REAL satu(ngrid,nlayer) ! satu ratio for output
	484	REAL zqsat(ngrid,nlayer) ! saturation
[42]	485
[414]	486	! Added for new NLTE scheme
[1047]	487	real co2vmr_gcm(ngrid,nlayer)
	488	real n2vmr_gcm(ngrid,nlayer)
	489	real ovmr_gcm(ngrid,nlayer)
	490	real covmr_gcm(ngrid,nlayer)
[1124]	491	integer ierr_nlte
	492	real*8 varerr
[414]	493
[2149]	494	c Non-oro GW tendencies
	495	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
	496	REAL d_t_hin(ngrid,nlayer)
[3144]	497	REAL d_u_mix(ngrid,nlayer), d_v_mix(ngrid,nlayer)
	498	REAL d_t_mix(ngrid,nlayer), zdq_mix(ngrid,nlayer,nq)
	499
[2149]	500	c Diagnostics 2D of gw_nonoro
	501	REAL zustrhi(ngrid), zvstrhi(ngrid)
[267]	502	c Variables for PBL
[1047]	503	REAL zz1(ngrid)
[1236]	504	REAL lmax_th_out(ngrid)
[1047]	505	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
	506	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
[3203]	507	INTEGER lmax_th(ngrid),n_out,n
[566]	508	CHARACTER(50) zstring
[1047]	509	REAL dtke_th(ngrid,nlayer+1)
[636]	510	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
[566]	511	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
[1047]	512	REAL u_out1(ngrid)
	513	REAL T_out1(ngrid)
[566]	514	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
[1236]	515	REAL tstar(ngrid) ! friction velocity and friction potential temp
[3203]	516	REAL vhf(ngrid), vvv(ngrid)
[1974]	517	real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer)
[3203]	518	real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)
[1974]	519	real qdusttotal0(ngrid),qdusttotal1(ngrid)
[790]	520
[1711]	521	c sub-grid scale water ice clouds (A. Pottier 2013)
	522	logical clearsky
	523	! flux for the part without clouds
[1972]	524	real zdtswclf(ngrid,nlayer)
	525	real zdtlwclf(ngrid,nlayer)
[3203]	526	real fluxsurf_lwclf(ngrid)
[2685]	527	real fluxsurf_dn_swclf(ngrid,2),fluxsurf_up_swclf(ngrid,2)
[1972]	528	real fluxtop_lwclf(ngrid)
[2685]	529	real fluxtop_dn_swclf(ngrid,2),fluxtop_up_swclf(ngrid,2)
[1972]	530	real taucloudtesclf(ngrid)
[1973]	531	real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds
[1972]	532	real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
[2362]	533	C test de conservation de la masse de CO2
	534	REAL co2totA
	535	REAL co2totB
[2660]	536	REAL co2conservation
[1711]	537
[2628]	538	c entrainment by mountain top dust flows above sub-grid scale topography
[2199]	539	REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons
	540
[2285]	541	c when no startfi file is asked for init
	542	real alpha,lay1 ! coefficients for building layers
	543	integer iloop
	544
[2551]	545	! flags to trigger extra sanity checks
	546	logical,save :: check_physics_inputs=.false.
	547	logical,save :: check_physics_outputs=.false.
[2281]	548
[3203]	549	!$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs)
	550
[2896]	551	c Sub-grid scale slopes
[2900]	552	real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K]
	553	real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1]
	554	real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1]
	555	real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2]
	556	real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie
[2896]	557	integer :: islope
[2900]	558	real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
[2896]	559
[2646]	560	logical :: write_restart
[2616]	561
[2963]	562	! Variable for ice table
[3230]	563	REAL :: rhowater_surf(ngrid,nslope) ! Water density at the surface [kg/m^3]
	564	REAL :: rhowater_surf_sat(ngrid,nslope) ! Water density at the surface at saturation [kg/m^3]
	565	REAL :: rhowater_soil(ngrid,nsoilmx,nslope) ! Water density in soil layers [kg/m^3]
	566	REAL,PARAMETER :: alpha_clap_h2o = 28.9074 ! Coeff for Clapeyron law [/]
	567	REAL,PARAMETER :: beta_clap_h2o = -6143.7 ! Coeff for Clapeyron law [K]
	568	REAL :: pvap_surf(ngrid) ! Water vapor partial pressure in first layer [Pa]
	569	REAL,PARAMETER :: m_co2 = 44.01E-3 ! CO2 molecular mass [kg/mol]
	570	REAL,PARAMETER :: m_noco2 = 33.37E-3 ! Non condensible mol mass [kg/mol]
	571	REAL :: ztmp1,ztmp2 ! intermediate variables to compute the mean molar mass of the layer
	572	REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
[3167]	573	! Variable for the computation of the TKE with parameterization from ATKE working group
	574	REAL :: viscom ! kinematic molecular viscosity for momentum
	575	REAL :: viscoh ! kinematic molecular viscosity for heat
	576
[42]	577	c=======================================================================
[2362]	578	pdq(:,:,:) = 0.
[42]	579
	580	c 1. Initialisation:
	581	c -----------------
	582	c 1.1 Initialisation only at first call
	583	c ---------------------------------------
	584	IF (firstcall) THEN
	585
[2551]	586	call getin_p("check_physics_inputs",check_physics_inputs)
	587	call getin_p("check_physics_outputs",check_physics_outputs)
	588
[42]	589	c variables set to 0
	590	c ~~~~~~~~~~~~~~~~~~
[286]	591	aerosol(:,:,:)=0
	592	dtrad(:,:)=0
[674]	593
	594	#ifndef MESOSCALE
[2900]	595	fluxrad(:,:)=0
[528]	596	wstar(:)=0.
[674]	597	#endif
[268]	598
[2223]	599	#ifdef CPP_XIOS
	600	! Initialize XIOS context
	601	write(,) "physiq: call wxios_context_init"
	602	CALL wxios_context_init
	603	#endif
	604
[3203]	605	c read startfi
[42]	606	c ~~~~~~~~~~~~
[226]	607	#ifndef MESOSCALE
[2896]	608
[1233]	609	! GCM. Read netcdf initial physical parameters.
[226]	610	CALL phyetat0 ("startfi.nc",0,0,
[3113]	611	& nsoilmx,ngrid,nlayer,nq,nqsoil,
[226]	612	& day_ini,time_phys,
[1944]	613	& tsurf,tsoil,albedo,emis,
[3113]	614	& q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,
[3130]	615	& watercap,perennial_co2ice,
[2999]	616	& def_slope,def_slope_mean,subslope_dist)
[185]	617
[3203]	618	! Sky view:
[2896]	619	DO islope=1,nslope
	620	sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
	621	END DO
	622	! Determine the 'flatest' slopes
	623	iflat = 1
	624	DO islope=2,nslope
	625	IF(abs(def_slope_mean(islope)).lt.
	626	& abs(def_slope_mean(iflat)))THEN
	627	iflat = islope
	628	ENDIF
	629	ENDDO
[3037]	630	write(,)'Flat slope for islope = ',iflat
	631	write(,)'corresponding criterium = ',def_slope_mean(iflat)
[2896]	632
[1233]	633	#else
	634	! MESOSCALE. Supposedly everything is already set in modules.
[1579]	635	! So we just check. And we fill day_ini
[3037]	636	write(,)"check: --- in physiq.F"
	637	write(,)"check: rad,cpp,g,r,rcp,daysec"
	638	write(,)rad,cpp,g,r,rcp,daysec
	639	write(,)'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
	640	write(,)'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
[3203]	641	write(,)'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
[3037]	642	write(,)'check: midlayer,layer ', mlayer(:),layer(:)
	643	write(,)'check: tracernames ', noms
	644	write(,)'check: emis ',emis(1,:),emis(ngrid,:)
	645	write(,)'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
	646	write(,)'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:)
	647	write(,)'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:)
[1579]	648	!!!
[1233]	649	day_ini = pday
[2491]	650	!!! a couple initializations (dummy for mesoscale) done in phyetat0
	651	!!! --- maybe this should be done in update_inputs_physiq_mod
[2562]	652
[2491]	653	tauscaling(:)=1.0 !! probably important
	654	totcloudfrac(:)=1.0
[2942]	655	DO islope = 1,nslope
	656	albedo(:,1,islope)=albedodat(:)
	657	albedo(:,2,islope)=albedo(:,1,islope)
	658	inertiesoil(:,:,islope) = inertiedat(:,:)
	659	watercap(:,:)=0.0
	660	ENDDO
[1233]	661	#endif
[2281]	662	#ifndef MESOSCALE
	663	if (.not.startphy_file) then
	664	! starting without startfi.nc and with callsoil
	665	! is not yet possible as soildepth default is not defined
[2362]	666	if (callsoil) then
[2285]	667	! default mlayer distribution, following a power law:
	668	! mlayer(k)=lay1alpha*(k-1/2)
	669	lay1=2.e-4
[3203]	670	alpha=2
[2285]	671	do iloop=0,nsoilmx-1
[3203]	672	mlayer(iloop)=lay1(alpha*(iloop-0.5))
	673	enddo
[2285]	674	lay1=sqrt(mlayer(0)*mlayer(1))
	675	alpha=mlayer(1)/mlayer(0)
	676	do iloop=1,nsoilmx
	677	layer(iloop)=lay1(alpha*(iloop-1))
	678	enddo
[2281]	679	endif
	680	! additionnal "academic" initialization of physics
[2900]	681	do islope = 1,nslope
	682	tsurf(:,islope)=pt(:,1)
	683	enddo
[2281]	684	write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	685	do isoil=1,nsoilmx
[2900]	686	tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:)
[2281]	687	enddo
	688	write(,) "Physiq: initializing inertiedat !!"
	689	inertiedat(:,:)=400.
[2942]	690	inertiesoil(:,:,:)=400.
[2281]	691	write(,) "Physiq: initializing day_ini to pdat !"
	692	day_ini=pday
	693	endif
	694	#endif
	695	if (pday.ne.day_ini) then
	696	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	697	& "physics and dynamics"
	698	write(,) "dynamics day [pday]: ",pday
	699	write(,) "physics day [day_ini]: ",day_ini
	700	call abort_physic("physiq","dynamics day /= physics day",1)
	701	endif
	702
	703	write (,) 'In physiq day_ini =', day_ini
	704
[286]	705	c initialize tracers
	706	c ~~~~~~~~~~~~~~~~~~
[2823]	707	CALL initracer(ngrid,nq,qsurf)
[286]	708
[42]	709	c Initialize albedo and orbital calculation
	710	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3013]	711	CALL surfini(ngrid,nslope,qsurf)
[42]	712	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[3203]	713	c initialize soil
[42]	714	c ~~~~~~~~~~~~~~~
	715	IF (callsoil) THEN
[833]	716	c Thermal inertia feedback:
[3230]	717	IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
[2900]	718	DO islope = 1,nslope
[3230]	719	CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
	720	s qsurf(:,igcm_h2o_ice,:),pore_icefraction,
	721	s inertiesoil_tifeedback)
[2900]	722	ENDDO
[2942]	723	CALL soil(ngrid,nsoilmx,firstcall,
	724	s inertiesoil_tifeedback,
[833]	725	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	726	ELSE
[2942]	727	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
[833]	728	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	729	ENDIF ! of IF (tifeedback)
[42]	730	ELSE
[3037]	731	write(,)
[42]	732	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	733	DO ig=1,ngrid
[2900]	734	capcal(ig,:)=1.e5
	735	fluxgrd(ig,:)=0.
[42]	736	ENDDO
	737	ENDIF
	738	icount=1
	739
[226]	740	#ifndef MESOSCALE
	741	c Initialize thermospheric parameters
	742	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[42]	743
[705]	744	if (callthermos) then
[1266]	745	call fill_data_thermos
	746	call allocate_param_thermos(nlayer)
	747	call allocate_param_iono(nlayer,nreact)
[1684]	748	call param_read_e107
[705]	749	endif
[226]	750	#endif
[42]	751
[3185]	752	c Initialize rnew cpnew and mmean as constant
	753	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3203]	754	call init_r_cp_mu(ngrid,nlayer)
[42]	755
[2383]	756	if(callnlte.and.nltemodel.eq.2) call nlte_setup
	757	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
[414]	758
[2383]	759
[2823]	760	IF (water.AND.(ngrid.NE.1)) THEN
[3203]	761	write(,)"physiq: water_param Surface water frost albedo:",
[2508]	762	. albedo_h2o_frost
[3203]	763	write(,)"physiq: water_param Surface watercap albedo:",
[2508]	764	. albedo_h2o_cap
[42]	765	ENDIF
[900]	766
	767	#ifndef MESOSCALE
[3203]	768	! no need to compute slopes when in 1D; it is an input
	769	if (ngrid /= 1 .and. callslope) call getslopes(ngrid,phisfi)
	770	if (ecritstart.GT.0) then
[2507]	771	call physdem0("restartfi.nc",longitude,latitude,
[2224]	772	& nsoilmx,ngrid,nlayer,nq,
[2544]	773	& ptimestep,pday,0.,cell_area,
[2896]	774	& albedodat,inertiedat,def_slope,
	775	& subslope_dist)
[3203]	776	else
[2507]	777	call physdem0("restartfi.nc",longitude,latitude,
	778	& nsoilmx,ngrid,nlayer,nq,
[2544]	779	& ptimestep,float(day_end),0.,cell_area,
[2896]	780	& albedodat,inertiedat,def_slope,
	781	& subslope_dist)
[3203]	782	endif
[2199]	783
[2628]	784	c Initialize mountain mesh fraction for the entrainment by top flows param.
[2199]	785	c ~~~~~~~~~~~~~~~
[2628]	786	if (topflows) call topmons_setup(ngrid)
[3203]	787
[3167]	788	c Parameterization of the ATKE routine
	789	c ~~~~~~~~~~~~~~~
	790	if (callatke) then
	791	viscom = 0.001
	792	viscoh = 0.001
	793	CALL atke_ini(g, r, pi, cpp, 0., viscom, viscoh)
	794	endif
[2212]	795
	796	#endif
[3203]	797
[2223]	798	#ifdef CPP_XIOS
	799	! XIOS outputs
	800	write(,) "physiq firstcall: call initialize_xios_output"
	801	call initialize_xios_output(pday,ptime,ptimestep,daysec,
[2333]	802	& presnivs,pseudoalt,mlayer)
[2223]	803	#endif
[42]	804	ENDIF ! (end of "if firstcall")
	805
[2551]	806	if (check_physics_inputs) then
	807	! Check the validity of input fields coming from the dynamics
[2570]	808	call check_physics_fields("begin physiq:",pt,pu,pv,pplev,pq)
[2551]	809	endif
	810
[42]	811	c ---------------------------------------------------
	812	c 1.2 Initializations done at every physical timestep:
	813	c ---------------------------------------------------
	814	c
[3203]	815	#ifdef CPP_XIOS
[2223]	816	! update XIOS time/calendar
[3203]	817	call update_xios_timestep
	818	#endif
[2312]	819
	820	c Initialize various variables
[42]	821	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[286]	822	pdv(:,:)=0
	823	pdu(:,:)=0
	824	pdt(:,:)=0
	825	pdq(:,:,:)=0
	826	pdpsrf(:)=0
[2900]	827	zflubid(:,:)=0
	828	zdtsurf(:,:)=0
	829	dqsurf(:,:,:)=0
[2252]	830	dsodust(:,:)=0.
	831	dsords(:,:)=0.
	832	dsotop(:,:)=0.
[2900]	833	dwatercap(:,:)=0
[2907]	834
	835	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
	836	& albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
[3203]	837
[2643]	838	! Dust scenario conversion coefficient from IRabs to VISext
	839	IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015
	840	! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true.
[2362]	841
[1377]	842	#ifdef DUSTSTORM
[1375]	843	pq_tmp(:,:,:)=0
[1377]	844	#endif
[3203]	845	igout=ngrid/2+1
[42]	846
	847
	848	zday=pday+ptime ! compute time, in sols (and fraction thereof)
[1974]	849	! Compute local time at each grid point
	850	DO ig=1,ngrid
	851	local_time(ig)=modulo(1.+(zday-INT(zday))
	852	& +(longitude_deg(ig)/15)/24,1.)
	853	ENDDO
[42]	854	c Compute Solar Longitude (Ls) :
	855	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	856	if (season) then
	857	call solarlong(zday,zls)
	858	else
	859	call solarlong(float(day_ini),zls)
	860	end if
	861
[883]	862	c Initialize pressure levels
	863	c ~~~~~~~~~~~~~~~~~~
	864	zplev(:,:) = pplev(:,:)
	865	zplay(:,:) = pplay(:,:)
	866	ps(:) = pplev(:,1)
	867
[3157]	868	#ifndef MESOSCALE
	869	c-----------------------------------------------------------------------
	870	c 1.2.1 Compute mean mass, cp, and R
[3185]	871	c update_r_cp_mu_ak outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
[3203]	872	c , mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
[3157]	873	c --------------------------------
	874
	875	if(photochem.or.callthermos) then
[3185]	876	call update_r_cp_mu_ak(ngrid,nlayer,nq,
[3157]	877	& zplay,pt,pdt,pq,pdq,ptimestep)
	878	endif
	879	#endif
	880
[42]	881	c Compute geopotential at interlayers
	882	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	883	c ponderation des altitudes au niveau des couches en dp/p
[3157]	884	cc ------------------------------------------
[3203]	885	!Calculation zzlev & zzlay for first layer
[42]	886	DO ig=1,ngrid
[3203]	887	zzlay(ig,1)=-(log(pplay(ig,1)/ps(ig)))rnew(ig,1)pt(ig,1)/g
	888	zzlev(ig,1)=0
[42]	889	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
[3157]	890	gz(ig,1)=g
[3203]	891
	892	DO l=2,nlayer
[3157]	893	! compute "mean" temperature of the layer
	894	if(pt(ig,l) .eq. pt(ig,l-1)) then
	895	tlaymean=pt(ig,l)
	896	else
	897	tlaymean=(pt(ig,l)- pt(ig,l-1))/log(pt(ig,l)/pt(ig,l-1))
	898	endif
[3203]	899
[3157]	900	! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
	901	gz(ig,l)=g(rad2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))*2
[3203]	902
	903	zzlay(ig,l)=zzlay(ig,l-1)-
[3157]	904	& (log(pplay(ig,l)/pplay(ig,l-1))rnew(ig,l)tlaymean/gz(ig,l))
[3203]	905
[3157]	906	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
	907	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
	908	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	909	ENDDO
[42]	910	ENDDO
	911
	912	! Potential temperature calculation not the same in physiq and dynamic
	913
	914	c Compute potential temperature
	915	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	916	DO l=1,nlayer
[3203]	917	DO ig=1,ngrid
[883]	918	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
[42]	919	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	920	ENDDO
	921	ENDDO
[2312]	922
[42]	923
[1313]	924	! Compute vertical velocity (m/s) from vertical mass flux
[1346]	925	! w = F / (rhoarea) and rho = P/(rT)
[1313]	926	! but first linearly interpolate mass flux to mid-layers
	927	do l=1,nlayer-1
[1312]	928	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
[1313]	929	enddo
	930	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	931	do l=1,nlayer
[3203]	932	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) /
[1541]	933	& (pplay(1:ngrid,l)*cell_area(1:ngrid))
[1346]	934	! NB: here we use r and not rnew since this diagnostic comes
[1312]	935	! from the dynamics
[1313]	936	enddo
[1312]	937
[2374]	938	! test for co2 conservation with co2 microphysics
	939	if (igcm_co2_ice.ne.0) then
	940	! calculates the amount of co2 at the beginning of physics
	941	co2totA = 0.
	942	do ig=1,ngrid
	943	do l=1,nlayer
	944	co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
	945	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
	946	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
	947	end do
[2900]	948	do islope = 1,nslope
	949	co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
	950	& subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
	951	enddo
[2362]	952	end do
[2599]	953	else
	954	co2totA = 0.
	955	do ig=1,ngrid
	956	do l=1,nlayer
	957	co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
	958	& (pq(ig,l,igcm_co2)
	959	& +pdq(ig,l,igcm_co2)*ptimestep)
	960	end do
[2900]	961	do islope = 1,nslope
	962	co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
	963	& subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
	964	enddo
[2599]	965	end do
[2374]	966	endif ! of if (igcm_co2_ice.ne.0)
[42]	967	c-----------------------------------------------------------------------
	968	c 2. Compute radiative tendencies :
	969	c------------------------------------
	970
	971	IF (callrad) THEN
	972
[2162]	973	c Local Solar zenith angle
	974	c ~~~~~~~~~~~~~~~~~~~~~~~~
[2880]	975
[2162]	976	CALL orbite(zls,dist_sol,declin)
[2880]	977
[2162]	978	IF (diurnal) THEN
	979	ztim1=SIN(declin)
	980	ztim2=COS(declin)COS(2.pi*(zday-.5))
	981	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
[42]	982
[2162]	983	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	984	& ztim1,ztim2,ztim3, mu0,fract)
[42]	985
[2162]	986	ELSE
	987	CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
	988	ENDIF ! of IF (diurnal)
[2880]	989
[2162]	990	IF( MOD(icount-1,iradia).EQ.0) THEN
	991
[42]	992	c NLTE cooling from CO2 emission
	993	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[414]	994	IF(callnlte) then
	995	if(nltemodel.eq.0.or.nltemodel.eq.1) then
[2447]	996	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
[414]	997	else if(nltemodel.eq.2) then
[705]	998	co2vmr_gcm(1:ngrid,1:nlayer)=
	999	& pq(1:ngrid,1:nlayer,igcm_co2)*
	1000	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
	1001	n2vmr_gcm(1:ngrid,1:nlayer)=
	1002	& pq(1:ngrid,1:nlayer,igcm_n2)*
	1003	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
	1004	covmr_gcm(1:ngrid,1:nlayer)=
	1005	& pq(1:ngrid,1:nlayer,igcm_co)*
	1006	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
	1007	ovmr_gcm(1:ngrid,1:nlayer)=
	1008	& pq(1:ngrid,1:nlayer,igcm_o)*
	1009	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
[2362]	1010
[883]	1011	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
[3203]	1012	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
[1124]	1013	$ ovmr_gcm, zdtnlte,ierr_nlte,varerr )
	1014	if(ierr_nlte.gt.0) then
	1015	write(,)
	1016	$ 'WARNING: nlte_tcool output with error message',
	1017	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
	1018	write(,)'I will continue anyway'
	1019	endif
[42]	1020
[705]	1021	zdtnlte(1:ngrid,1:nlayer)=
	1022	& zdtnlte(1:ngrid,1:nlayer)/86400.
[414]	1023	endif
[1974]	1024	ELSE
[528]	1025	zdtnlte(:,:)=0.
[1974]	1026	ENDIF !end callnlte
[42]	1027
[3369]	1028	! Find number of layers for LTE radiation calculations
	1029	! (done only once per sol)
	1030	IF(MOD((icount-1),steps_per_sol).EQ.0)
[883]	1031	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
[42]	1032
[1974]	1033	c rocketstorm : compute dust storm mesh fraction
	1034	IF (rdstorm) THEN
	1035	CALL calcstormfract(ngrid,nlayer,nq,pq,
	1036	& totstormfract)
	1037	ENDIF
	1038
[42]	1039	c Note: Dustopacity.F has been transferred to callradite.F
[1410]	1040
	1041	#ifdef DUSTSTORM
	1042	!! specific case: save the quantity of dust before adding perturbation
[2616]	1043
[1410]	1044	if (firstcall) then
	1045	pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
	1046	pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
	1047	endif
	1048	#endif
[3203]	1049
[42]	1050	c Call main radiative transfer scheme
	1051	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1052	c Transfer through CO2 (except NIR CO2 absorption)
	1053	c and aerosols (dust and water ice)
[2628]	1054	! callradite for background dust (out of the rdstorm fraction)
[1974]	1055	clearatm=.true.
[2628]	1056	!! callradite for background dust (out of the topflows fraction)
[2199]	1057	nohmons=.true.
[2312]	1058
[42]	1059	c Radiative transfer
	1060	c ------------------
[1711]	1061	! callradite for the part with clouds
[2362]	1062	clearsky=.false. ! part with clouds for both cases CLFvarying true/false
[2900]	1063	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
[2907]	1064	& albedo_meshavg,emis_meshavg,
[2908]	1065	& mu0,zplev,zplay,pt,tsurf_meshavg,fract,dist_sol,igout,
[2900]	1066	& zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
	1067	& fluxsurf_up_sw,
[2685]	1068	& fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
	1069	& tau_pref_scenario,tau_pref_gcm,
[2643]	1070	& tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
[2826]	1071	& taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
[2907]	1072	& qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract,
[2826]	1073	& clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
[2616]	1074
[2900]	1075	DO islope=1,nslope
	1076	fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat)
	1077	fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat)
	1078	ENDDO
	1079
[1711]	1080	! case of sub-grid water ice clouds: callradite for the clear case
	1081	IF (CLFvarying) THEN
[3203]	1082	! ---> PROBLEMS WITH ALLOCATED ARRAYS
[1711]	1083	! (temporary solution in callcorrk: do not deallocate
	1084	! if
	1085	! CLFvarying ...) ?? AP ??
[3203]	1086	clearsky=.true.
[1711]	1087	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
[2907]	1088	& albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt,
[2908]	1089	& tsurf_meshavg,fract,
[2685]	1090	& dist_sol,igout,zdtlwclf,zdtswclf,
	1091	& fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf,
	1092	& fluxtop_lwclf,fluxtop_dn_swclf,fluxtop_up_swclf,
[2643]	1093	& tau_pref_scenario,tau_pref_gcm,tau,aerosol,
	1094	& dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
[2417]	1095	& taucloudtesclf,rdust,
[2907]	1096	& rice,nuice,riceco2, nuiceco2,
	1097	& qsurf_meshavg(:,igcm_co2),
[2826]	1098	& rstormdust,rtopdust,totstormfract,
[2628]	1099	& clearatm,dsords,dsotop,
	1100	& nohmons,clearsky,totcloudfrac)
[1711]	1101	clearsky = .false. ! just in case.
	1102	! Sum the fluxes and heating rates from cloudy/clear
	1103	! cases
	1104	DO ig=1,ngrid
[1973]	1105	tf_clf=totcloudfrac(ig)
	1106	ntf_clf=1.-tf_clf
[2900]	1107	DO islope=1,nslope
[3203]	1108	fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig)
[2900]	1109	& + tf_clf*fluxsurf_lw(ig,islope)
[3203]	1110	fluxsurf_dn_sw(ig,1:2,islope) =
	1111	& ntf_clf*fluxsurf_dn_swclf(ig,1:2)
[2900]	1112	& + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
	1113	ENDDO
[3203]	1114	fluxsurf_up_sw(ig,1:2) =
	1115	& ntf_clf*fluxsurf_up_swclf(ig,1:2)
[2685]	1116	& + tf_clf*fluxsurf_up_sw(ig,1:2)
[3203]	1117	fluxtop_lw(ig) = ntf_clf*fluxtop_lwclf(ig)
[1973]	1118	& + tf_clf*fluxtop_lw(ig)
[3203]	1119	fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2)
[2685]	1120	& + tf_clf*fluxtop_dn_sw(ig,1:2)
[3203]	1121	fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2)
[2685]	1122	& + tf_clf*fluxtop_up_sw(ig,1:2)
[3203]	1123	taucloudtes(ig) = ntf_clf*taucloudtesclf(ig)
[1973]	1124	& + tf_clf*taucloudtes(ig)
[3203]	1125	zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer)
[1973]	1126	& + tf_clf*zdtlw(ig,1:nlayer)
[3203]	1127	zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer)
[1973]	1128	& + tf_clf*zdtsw(ig,1:nlayer)
[1711]	1129	ENDDO
[42]	1130
[1711]	1131	ENDIF ! (CLFvarying)
[3203]	1132
[2265]	1133	!============================================================================
[3203]	1134
[1375]	1135	#ifdef DUSTSTORM
[1410]	1136	!! specific case: compute the added quantity of dust for perturbation
	1137	if (firstcall) then
[3203]	1138	pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
	1139	& pdq(1:ngrid,1:nlayer,igcm_dust_mass)
[1410]	1140	& - pq_tmp(1:ngrid,1:nlayer,1)
	1141	& + pq(1:ngrid,1:nlayer,igcm_dust_mass)
	1142	pdq(1:ngrid,1:nlayer,igcm_dust_number)=
	1143	& pdq(1:ngrid,1:nlayer,igcm_dust_number)
	1144	& - pq_tmp(1:ngrid,1:nlayer,2)
	1145	& + pq(1:ngrid,1:nlayer,igcm_dust_number)
	1146	endif
[1375]	1147	#endif
	1148
[234]	1149	c Outputs for basic check (middle of domain)
	1150	c ------------------------------------------
[627]	1151	write(*,'("Ls =",f11.6," check lat =",f10.6,
	1152	& " lon =",f11.6)')
[1541]	1153	& zls180./pi,latitude(igout)180/pi,
	1154	& longitude(igout)*180/pi
[2281]	1155
[2415]	1156	write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6,
[627]	1157	& " tau(",f4.0," Pa) =",f9.6)')
[2415]	1158	& odpref,tau_pref_gcm(igout),
[883]	1159	& odpref,tau(igout,1)*odpref/zplev(igout,1)
[2616]	1160
	1161
[234]	1162	c ---------------------------------------------------------
	1163	c Call slope parameterization for direct and scattered flux
	1164	c ---------------------------------------------------------
	1165	IF(callslope) THEN
[2900]	1166	! assume that in this case, nslope = 1
	1167	if(nslope.ne.1) then
	1168	call abort_physic(
	1169	& "physiq","callslope=true but nslope.ne.1",1)
	1170	endif
[3203]	1171	write(,) 'Slope scheme is on and computing...'
	1172	DO ig=1,ngrid
[234]	1173	sl_the = theta_sl(ig)
	1174	IF (sl_the .ne. 0.) THEN
[2942]	1175	ztim1=fluxsurf_dn_sw(ig,1,iflat)
	1176	& +fluxsurf_dn_sw(ig,2,iflat)
[234]	1177	DO l=1,2
[1541]	1178	sl_lct = ptime24. + 180.longitude(ig)/pi/15.
[234]	1179	sl_ra = pi*(1.0-sl_lct/12.)
[1541]	1180	sl_lat = 180.*latitude(ig)/pi
[577]	1181	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
[2942]	1182	sl_alb = albedo(ig,l,iflat)
[234]	1183	sl_psi = psi_sl(ig)
[2942]	1184	sl_fl0 = fluxsurf_dn_sw(ig,l,iflat)
[234]	1185	sl_di0 = 0.
[2879]	1186	if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
[3203]	1187	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
[2678]	1188	sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
[234]	1189	sl_di0 = sl_di0/ztim1
[2942]	1190	sl_di0 = fluxsurf_dn_sw(ig,l,iflat)*sl_di0
[234]	1191	endif
	1192	! you never know (roundup concern...)
	1193	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
	1194	!!!!!!!!!!!!!!!!!!!!!!!!!!
[3203]	1195	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
[234]	1196	& sl_tau, sl_alb, sl_the, sl_psi,
	1197	& sl_di0, sl_fl0, sl_flu )
	1198	!!!!!!!!!!!!!!!!!!!!!!!!!!
[2900]	1199	fluxsurf_dn_sw(ig,l,1) = sl_flu
[234]	1200	ENDDO
	1201	!!! compute correction on IR flux as well
	1202	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
[2900]	1203	fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky
[234]	1204	ENDIF
	1205	ENDDO
[2942]	1206	ELSE ! not calling subslope, nslope might be > 1
	1207	DO islope = 1,nslope
[3203]	1208	sl_the=abs(def_slope_mean(islope))
[2942]	1209	IF (sl_the .gt. 1e-6) THEN
[2953]	1210	IF(def_slope_mean(islope).ge.0.) THEN
	1211	psi_sl(:) = 0. !Northward slope
	1212	ELSE
	1213	psi_sl(:) = 180. !Southward slope
	1214	ENDIF
[3203]	1215	DO ig=1,ngrid
[2942]	1216	ztim1=fluxsurf_dn_sw(ig,1,islope)
	1217	s +fluxsurf_dn_sw(ig,2,islope)
	1218	DO l=1,2
	1219	sl_lct = ptime24. + 180.longitude(ig)/pi/15.
	1220	sl_ra = pi*(1.0-sl_lct/12.)
	1221	sl_lat = 180.*latitude(ig)/pi
	1222	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
	1223	sl_alb = albedo(ig,l,islope)
	1224	sl_psi = psi_sl(ig)
	1225	sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
	1226	sl_di0 = 0.
	1227	if (mu0(ig) .gt. 0.) then
[3203]	1228	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
[2942]	1229	sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
	1230	sl_di0 = sl_di0/ztim1
	1231	sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
	1232	endif
	1233	! you never know (roundup concern...)
	1234	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
	1235	!!!!!!!!!!!!!!!!!!!!!!!!!!
[3203]	1236	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
[2942]	1237	& sl_tau, sl_alb, sl_the, sl_psi,
	1238	& sl_di0, sl_fl0, sl_flu )
	1239	!!!!!!!!!!!!!!!!!!!!!!!!!!
	1240	fluxsurf_dn_sw(ig,l,islope) = sl_flu
	1241	ENDDO
	1242	!!! compute correction on IR flux as well
	1243
	1244	fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
	1245	& *sky_slope(islope)
	1246	ENDDO
	1247	ENDIF ! sub grid is not flat
	1248	ENDDO ! islope = 1,nslope
[2900]	1249	ENDIF ! callslope
[234]	1250
[42]	1251	c CO2 near infrared absorption
	1252	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[286]	1253	zdtnirco2(:,:)=0
[42]	1254	if (callnirco2) then
[883]	1255	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
[42]	1256	. mu0,fract,declin, zdtnirco2)
	1257	endif
	1258
	1259	c Radiative flux from the sky absorbed by the surface (W.m-2)
	1260	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1261	DO ig=1,ngrid
[2900]	1262	DO islope = 1,nslope
	1263	fluxrad_sky(ig,islope) =
	1264	$ emis(ig,islope)*fluxsurf_lw(ig,islope)
[3203]	1265	$ +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))
[2900]	1266	$ +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
	1267	ENDDO
[42]	1268	ENDDO
	1269
	1270	c Net atmospheric radiative heating rate (K.s-1)
	1271	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1272	IF(callnlte) THEN
[883]	1273	CALL blendrad(ngrid, nlayer, zplay,
[42]	1274	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
	1275	ELSE
	1276	DO l=1,nlayer
	1277	DO ig=1,ngrid
	1278	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	1279	& +zdtnirco2(ig,l)
	1280	ENDDO
	1281	ENDDO
	1282	ENDIF
	1283
	1284	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	1285
	1286	c Transformation of the radiative tendencies:
	1287	c -------------------------------------------
	1288
	1289	c Net radiative surface flux (W.m-2)
	1290	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
[2953]	1291
[42]	1292	c
	1293	DO ig=1,ngrid
[2900]	1294	DO islope = 1,nslope
	1295	zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope)
	1296	zplanck(ig)=emis(ig,islope)*
[42]	1297	$ stephanzplanck(ig)zplanck(ig)
[2900]	1298	fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig)
[234]	1299	IF(callslope) THEN
	1300	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
[2900]	1301	fluxrad(ig,nslope)=fluxrad(ig,nslope)+
	1302	$ (1.-sky)*zplanck(ig)
[3203]	1303	ELSE
[2900]	1304	fluxrad(ig,islope)=fluxrad(ig,islope) +
	1305	$ (1.-sky_slope(iflat))emis(ig,iflat)
[3203]	1306	$ stephantsurf(ig,iflat)*4
[234]	1307	ENDIF
[42]	1308	ENDDO
[2900]	1309	ENDDO
[42]	1310
	1311	DO l=1,nlayer
	1312	DO ig=1,ngrid
	1313	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	1314	ENDDO
	1315	ENDDO
	1316
	1317	ENDIF ! of IF (callrad)
	1318
[2199]	1319	c 3.1 Rocket dust storm
[1974]	1320	c -------------------------------------------
	1321	IF (rdstorm) THEN
	1322	clearatm=.false.
	1323	pdqrds(:,:,:)=0.
	1324	qdusttotal0(:)=0.
	1325	qdusttotal1(:)=0.
	1326	do ig=1,ngrid
	1327	do l=1,nlayer
	1328	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential
	1329	! temperature tendency
	1330	! for diagnostics
[2934]	1331	! qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+
	1332	! & pdq(ig,l,igcm_dust_mass)*ptimestep
	1333	! qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+
	1334	! & pdq(ig,l,igcm_stormdust_mass)*ptimestep
	1335	! qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+
	1336	! & qstormrds0(ig,l))*(zplev(ig,l)-
	1337	! & zplev(ig,l+1))/g
[1974]	1338	enddo
	1339	enddo
[2934]	1340	! call write_output('qdustrds0','qdust before rds',
	1341	! & 'kg/kg ',qdustrds0(:,:))
	1342	! call write_output('qstormrds0','qstorm before rds',
	1343	! & 'kg/kg ',qstormrds0(:,:))
[1974]	1344
	1345	CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,
	1346	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
	1347	& zzlay,zdtsw,zdtlw,
	1348	c for radiative transfer
	1349	& clearatm,icount,zday,zls,
[2953]	1350	& tsurf_meshavg,qsurf_meshavg(:,igcm_co2),
	1351	& igout,totstormfract,tauscaling,
[2826]	1352	& dust_rad_adjust,IRtoVIScoef,
[2953]	1353	& albedo_meshavg,emis_meshavg,
[1974]	1354	c input sub-grid scale cloud
	1355	& clearsky,totcloudfrac,
[2199]	1356	c input sub-grid scale topography
[2628]	1357	& nohmons,
[1974]	1358	c output
[2246]	1359	& pdqrds,wspeed,dsodust,dsords,dsotop,
[2415]	1360	& tau_pref_scenario,tau_pref_gcm)
[3203]	1361
[1974]	1362	c update the tendencies of both dust after vertical transport
	1363	DO l=1,nlayer
	1364	DO ig=1,ngrid
	1365	pdq(ig,l,igcm_stormdust_mass)=
[3203]	1366	& pdq(ig,l,igcm_stormdust_mass)+
[1974]	1367	& pdqrds(ig,l,igcm_stormdust_mass)
	1368	pdq(ig,l,igcm_stormdust_number)=
	1369	& pdq(ig,l,igcm_stormdust_number)+
	1370	& pdqrds(ig,l,igcm_stormdust_number)
	1371
	1372	pdq(ig,l,igcm_dust_mass)=
	1373	& pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
	1374	pdq(ig,l,igcm_dust_number)=
[3203]	1375	& pdq(ig,l,igcm_dust_number)+
[1974]	1376	& pdqrds(ig,l,igcm_dust_number)
	1377
	1378	ENDDO
	1379	ENDDO
	1380	do l=1,nlayer
	1381	do ig=1,ngrid
	1382	qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
	1383	& pdq(ig,l,igcm_dust_mass)*ptimestep
	1384	qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
	1385	& pdq(ig,l,igcm_stormdust_mass)*ptimestep
	1386	qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
	1387	& qstormrds1(ig,l))*(zplev(ig,l)-
	1388	& zplev(ig,l+1))/g
	1389	enddo
	1390	enddo
	1391
	1392	c for diagnostics
[2934]	1393	! call write_output('qdustrds1','qdust after rds',
	1394	! & 'kg/kg ',qdustrds1(:,:))
	1395	! call write_output('qstormrds1','qstorm after rds',
	1396	! & 'kg/kg ',qstormrds1(:,:))
[3203]	1397	!
[2934]	1398	! call write_output('qdusttotal0','q sum before rds',
	1399	! & 'kg/m2 ',qdusttotal0(:))
	1400	! call write_output('qdusttotal1','q sum after rds',
	1401	! & 'kg/m2 ',qdusttotal1(:))
[1974]	1402
	1403	ENDIF ! end of if(rdstorm)
	1404
[2199]	1405	c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography
	1406	c -------------------------------------------
[2628]	1407	IF (topflows) THEN
[2362]	1408	clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
[2199]	1409	nohmons=.false.
[2362]	1410	pdqtop(:,:,:)=0.
[2199]	1411	CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
	1412	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
	1413	& zzlay,zdtsw,zdtlw,
[2900]	1414	& icount,zday,zls,tsurf(:,iflat),
[3203]	1415	& qsurf_meshavg(:,igcm_co2),
[2826]	1416	& igout,aerosol,tauscaling,dust_rad_adjust,
[2953]	1417	& IRtoVIScoef,albedo_meshavg,emis_meshavg,
	1418	& totstormfract,clearatm,
[2199]	1419	& clearsky,totcloudfrac,
[2628]	1420	& nohmons,
[2246]	1421	& pdqtop,wtop,dsodust,dsords,dsotop,
[2415]	1422	& tau_pref_scenario,tau_pref_gcm)
[3203]	1423
[2199]	1424	c update the tendencies of both dust after vertical transport
	1425	DO l=1,nlayer
	1426	DO ig=1,ngrid
	1427	pdq(ig,l,igcm_topdust_mass)=
	1428	& pdq(ig,l,igcm_topdust_mass)+
	1429	& pdqtop(ig,l,igcm_topdust_mass)
	1430	pdq(ig,l,igcm_topdust_number)=
	1431	& pdq(ig,l,igcm_topdust_number)+
	1432	& pdqtop(ig,l,igcm_topdust_number)
	1433	pdq(ig,l,igcm_dust_mass)=
	1434	& pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass)
	1435	pdq(ig,l,igcm_dust_number)=
	1436	& pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number)
	1437
	1438	ENDDO
	1439	ENDDO
	1440
[2628]	1441	ENDIF ! end of if (topflows)
[2199]	1442
	1443	c 3.3 Dust injection from the surface
	1444	c -------------------------------------------
	1445	if (dustinjection.gt.0) then
[2281]	1446
[2199]	1447	CALL compute_dtau(ngrid,nlayer,
[3203]	1448	& zday,pplev,tau_pref_gcm,
	1449	& ptimestep,local_time,IRtoVIScoef,
[2643]	1450	& dustliftday)
[2199]	1451	endif ! end of if (dustinjection.gt.0)
	1452
[226]	1453	c-----------------------------------------------------------------------
[1974]	1454	c 4. Gravity wave and subgrid scale topography drag :
[226]	1455	c -------------------------------------------------
[42]	1456
[226]	1457
	1458	IF(calllott)THEN
	1459	CALL calldrag_noro(ngrid,nlayer,ptimestep,
[883]	1460	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
[2362]	1461
[226]	1462	DO l=1,nlayer
	1463	DO ig=1,ngrid
	1464	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
	1465	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
	1466	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
	1467	ENDDO
	1468	ENDDO
	1469	ENDIF
[234]	1470
[42]	1471	c-----------------------------------------------------------------------
[1974]	1472	c 5. Vertical diffusion (turbulent mixing):
[42]	1473	c -----------------------------------------
[226]	1474
[42]	1475	IF (calldifv) THEN
	1476	DO ig=1,ngrid
[2900]	1477	DO islope = 1,nslope
	1478	zflubid(ig,islope)=fluxrad(ig,islope)
	1479	& +fluxgrd(ig,islope)
	1480	ENDDO
[42]	1481	ENDDO
[286]	1482	zdum1(:,:)=0
	1483	zdum2(:,:)=0
[42]	1484	do l=1,nlayer
	1485	do ig=1,ngrid
	1486	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	1487	enddo
	1488	enddo
[226]	1489
[288]	1490	c ----------------------
[284]	1491	c Treatment of a special case : using new surface layer (Richardson based)
	1492	c without using the thermals in gcm and mesoscale can yield problems in
	1493	c weakly unstable situations when winds are near to 0. For those cases, we add
[3203]	1494	c a unit subgrid gustiness. Remember that thermals should be used we using the
[284]	1495	c Richardson based surface layer model.
[3203]	1496	IF ( .not.calltherm
	1497	. .and. callrichsl
[1236]	1498	. .and. .not.turb_resolved) THEN
[2616]	1499
[1047]	1500	DO ig=1, ngrid
[2908]	1501	IF (zh(ig,1) .lt. tsurf_meshavg(ig)) THEN
[528]	1502	wstar(ig)=1.
	1503	hfmax_th(ig)=0.2
	1504	ELSE
	1505	wstar(ig)=0.
	1506	hfmax_th(ig)=0.
[3203]	1507	ENDIF
[284]	1508	ENDDO
	1509	ENDIF
[2953]	1510
[288]	1511	c ----------------------
[284]	1512
[544]	1513	IF (tke_heat_flux .ne. 0.) THEN
[2616]	1514
[544]	1515	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
[883]	1516	& (-alog(zplay(:,1)/zplev(:,1)))
[544]	1517	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
	1518	ENDIF
[2312]	1519
[42]	1520	c Calling vdif (Martian version WITH CO2 condensation)
[2900]	1521	dwatercap_dif(:,:) = 0.
[3115]	1522	CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
[3468]	1523	$ ptimestep,capcal,
[883]	1524	$ zplay,zplev,zzlay,zzlev,z0,
[3230]	1525	$ pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,
	1526	$ qsoil,pore_icefraction,
[42]	1527	$ zdum1,zdum2,zdh,pdq,zflubid,
	1528	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
[3203]	1529	& zdqdif,zdqsdif,wstar,hfmax_th,
[1996]	1530	& zcondicea_co2microp,sensibFlux,
[3115]	1531	& dustliftday,local_time,watercap,dwatercap_dif)
[2896]	1532
[1236]	1533	DO ig=1,ngrid
[2900]	1534	zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:)
	1535	dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
[1236]	1536	ENDDO
[2963]	1537
[2907]	1538	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
	1539	& albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
[1236]	1540	IF (.not.turb_resolved) THEN
	1541	DO l=1,nlayer
[42]	1542	DO ig=1,ngrid
	1543	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	1544	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	1545	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	1546
	1547	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	1548	ENDDO
[226]	1549	ENDDO
[2312]	1550
[42]	1551	DO iq=1, nq
	1552	DO l=1,nlayer
	1553	DO ig=1,ngrid
[2362]	1554	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
[42]	1555	ENDDO
	1556	ENDDO
	1557	ENDDO
	1558	DO iq=1, nq
	1559	DO ig=1,ngrid
[2900]	1560	dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
[42]	1561	ENDDO
	1562	ENDDO
[2823]	1563
[1242]	1564	ELSE
	1565	write (,) '******************************************'
	1566	write (,) '** LES mode: the difv part is only used to'
	1567	write (,) '** - provide HFX and UST to the dynamics'
	1568	write (,) '** - update TSURF'
	1569	write (,) '******************************************'
	1570	!! Specific treatment for lifting in turbulent-resolving mode (AC)
	1571	IF (lifting .and. doubleq) THEN
[3203]	1572	!! lifted dust is injected in the first layer.
	1573	!! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
[1242]	1574	!! => lifted dust is not incremented before the sedimentation step.
	1575	zdqdif(1:ngrid,1,1:nq)=0.
[3203]	1576	zdqdif(1:ngrid,1,igcm_dust_number) =
[2907]	1577	. -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number)
[3203]	1578	zdqdif(1:ngrid,1,igcm_dust_mass) =
[2907]	1579	. -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass)
[1266]	1580	zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
[1242]	1581	DO iq=1, nq
[1974]	1582	IF ((iq .ne. igcm_dust_mass)
	1583	& .and. (iq .ne. igcm_dust_number)) THEN
[2900]	1584	zdqsdif(:,iq,:)=0.
[1974]	1585	ENDIF
[1242]	1586	ENDDO
	1587	ELSE
	1588	zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
[2900]	1589	zdqsdif(1:ngrid,1:nq,1:nslope) = 0.
[1242]	1590	ENDIF
[790]	1591	ENDIF
[3203]	1592	ELSE
[42]	1593	DO ig=1,ngrid
[2900]	1594	DO islope=1,nslope
	1595	zdtsurf(ig,islope)=zdtsurf(ig,islope)+
	1596	s (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope)
	1597	ENDDO
[42]	1598	ENDDO
[2963]	1599
[1236]	1600	IF (turb_resolved) THEN
[3203]	1601	write(,) 'Turbulent-resolving mode !'
[42]	1602	write(,) 'Please set calldifv to T in callphys.def'
[2214]	1603	call abort_physic("physiq","turbulent-resolving mode",1)
[42]	1604	ENDIF
	1605	ENDIF ! of IF (calldifv)
	1606
[162]	1607	c-----------------------------------------------------------------------
[1974]	1608	c 6. Thermals :
[162]	1609	c -----------------------------
[566]	1610
[1236]	1611	if(calltherm .and. .not.turb_resolved) then
[2362]	1612
[2823]	1613	call calltherm_interface(ngrid,nlayer,nq,igcm_co2,
[652]	1614	$ zzlev,zzlay,
[162]	1615	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
[883]	1616	$ zplay,zplev,pphi,zpopsk,
[185]	1617	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
[660]	1618	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
[3203]	1619
[162]	1620	DO l=1,nlayer
	1621	DO ig=1,ngrid
	1622	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
	1623	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
	1624	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
	1625	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
	1626	ENDDO
	1627	ENDDO
[2362]	1628
[162]	1629	DO ig=1,ngrid
	1630	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
[2362]	1631	ENDDO
[2312]	1632
[162]	1633	DO iq=1,nq
	1634	DO l=1,nlayer
	1635	DO ig=1,ngrid
	1636	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
	1637	ENDDO
	1638	ENDDO
	1639	ENDDO
[42]	1640
[277]	1641	lmax_th_out(:)=real(lmax_th(:))
	1642
[1236]	1643	else !of if calltherm
[162]	1644	lmax_th(:)=0
[528]	1645	wstar(:)=0.
	1646	hfmax_th(:)=0.
[277]	1647	lmax_th_out(:)=0.
[1236]	1648	end if
[2312]	1649
[42]	1650	c-----------------------------------------------------------------------
[1974]	1651	c 7. Dry convective adjustment:
[42]	1652	c -----------------------------
	1653
	1654	IF(calladj) THEN
	1655
	1656	DO l=1,nlayer
	1657	DO ig=1,ngrid
	1658	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	1659	ENDDO
	1660	ENDDO
[286]	1661	zduadj(:,:)=0
	1662	zdvadj(:,:)=0
	1663	zdhadj(:,:)=0
[42]	1664
	1665	CALL convadj(ngrid,nlayer,nq,ptimestep,
[883]	1666	$ zplay,zplev,zpopsk,lmax_th,
[42]	1667	$ pu,pv,zh,pq,
	1668	$ pdu,pdv,zdh,pdq,
	1669	$ zduadj,zdvadj,zdhadj,
	1670	$ zdqadj)
	1671
	1672	DO l=1,nlayer
	1673	DO ig=1,ngrid
	1674	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	1675	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	1676	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	1677
	1678	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	1679	ENDDO
	1680	ENDDO
	1681
	1682	DO iq=1, nq
	1683	DO l=1,nlayer
	1684	DO ig=1,ngrid
[3203]	1685	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
[42]	1686	ENDDO
	1687	ENDDO
	1688	ENDDO
	1689	ENDIF ! of IF(calladj)
	1690
[2149]	1691	c-----------------------------------------------------
	1692	c 8. Non orographic Gravity waves :
	1693	c -------------------------------------------------
[42]	1694
[2149]	1695	IF (calllott_nonoro) THEN
	1696
[2594]	1697	CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,
	1698	& cpnew,rnew,
	1699	& zplay,
[2149]	1700	& zmax_th, ! max altitude reached by thermals (m)
	1701	& pt, pu, pv,
	1702	& pdt, pdu, pdv,
	1703	& zustrhi,zvstrhi,
[2225]	1704	& d_t_hin, d_u_hin, d_v_hin)
[3144]	1705	IF (calljliu_gwimix) THEN
	1706	CALL nonoro_gwd_mix(ngrid,nlayer,ptimestep,
[3203]	1707	& nq,cpnew, rnew,
	1708	& zplay,
	1709	& zmax_th,
[3393]	1710	& pt, pu, pv, pq, zh,
[3144]	1711	!loss, chemical reaction loss rates
[3393]	1712	& pdt, pdu, pdv, pdq, zdh,
[3144]	1713	! zustrhi,zvstrhi,
	1714	& zdq_mix, d_t_mix, d_u_mix, d_v_mix)
	1715	ENDIF
[2149]	1716
	1717	! Update tendencies
	1718	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
	1719	& +d_t_hin(1:ngrid,1:nlayer)
	1720	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
[3203]	1721	& +d_u_hin(1:ngrid,1:nlayer)
[2149]	1722	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
	1723	& +d_v_hin(1:ngrid,1:nlayer)
[3144]	1724	! Update tendencies of gw mixing
	1725	IF (calljliu_gwimix) THEN
	1726	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
	1727	& +d_t_mix(1:ngrid,1:nlayer)
	1728	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
[3203]	1729	& +d_u_mix(1:ngrid,1:nlayer)
[3144]	1730	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
	1731	& +d_v_mix(1:ngrid,1:nlayer)
	1732	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)
	1733	& +zdq_mix(1:ngrid,1:nlayer,1:nq)
	1734	ENDIF
[2149]	1735
[3144]	1736
[2149]	1737	ENDIF ! of IF (calllott_nonoro)
	1738
[42]	1739	c-----------------------------------------------------------------------
[3203]	1740	c 9. Specific parameterizations for tracers
[42]	1741	c: -----------------------------------------
	1742
	1743
[2149]	1744	c 9a. Water and ice
[42]	1745	c ---------------
	1746
	1747	c ---------------------------------------
	1748	c Water ice condensation in the atmosphere
	1749	c ----------------------------------------
	1750	IF (water) THEN
	1751
	1752	call watercloud(ngrid,nlayer,ptimestep,
[883]	1753	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
[626]	1754	& pq,pdq,zdqcloud,zdtcloud,
[358]	1755	& nq,tau,tauscaling,rdust,rice,nuice,
[1711]	1756	& rsedcloud,rhocloud,totcloudfrac)
[633]	1757	c Temperature variation due to latent heat release
[42]	1758	if (activice) then
[2362]	1759	pdt(1:ngrid,1:nlayer) =
	1760	& pdt(1:ngrid,1:nlayer) +
[633]	1761	& zdtcloud(1:ngrid,1:nlayer)
[42]	1762	endif
[2312]	1763
[42]	1764	! increment water vapour and ice atmospheric tracers tendencies
[3203]	1765	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
	1766	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
[706]	1767	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
[3203]	1768	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
	1769	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
[706]	1770	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
	1771
[2312]	1772	if (hdo) then
	1773	! increment HDO vapour and ice atmospheric tracers tendencies
[2362]	1774	pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
	1775	& pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
[2312]	1776	& zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
[2362]	1777	pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
	1778	& pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
[2312]	1779	& zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
	1780	endif !hdo
	1781
[706]	1782	! increment dust and ccn masses and numbers
[883]	1783	! We need to check that we have Nccn & Ndust > 0
	1784	! This is due to single precision rounding problems
[706]	1785	if (microphys) then
[3203]	1786	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
	1787	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
[706]	1788	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
[3203]	1789	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
	1790	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
[706]	1791	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
[3203]	1792	where (pq(:,:,igcm_ccn_mass) +
[883]	1793	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
[3203]	1794	pdq(:,:,igcm_ccn_mass) =
[883]	1795	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
[3203]	1796	pdq(:,:,igcm_ccn_number) =
[883]	1797	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1798	end where
[3203]	1799	where (pq(:,:,igcm_ccn_number) +
[883]	1800	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
[3203]	1801	pdq(:,:,igcm_ccn_mass) =
[883]	1802	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
[3203]	1803	pdq(:,:,igcm_ccn_number) =
[883]	1804	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1805	end where
[706]	1806	endif
	1807
[883]	1808	if (scavenging) then
[3203]	1809	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
	1810	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
[883]	1811	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
[3203]	1812	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
	1813	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
[883]	1814	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
[3203]	1815	where (pq(:,:,igcm_dust_mass) +
[883]	1816	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
[3203]	1817	pdq(:,:,igcm_dust_mass) =
[883]	1818	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
[3203]	1819	pdq(:,:,igcm_dust_number) =
[883]	1820	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1821	end where
[3203]	1822	where (pq(:,:,igcm_dust_number) +
[883]	1823	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
[3203]	1824	pdq(:,:,igcm_dust_mass) =
[883]	1825	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
[3203]	1826	pdq(:,:,igcm_dust_number) =
[883]	1827	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1828	end where
	1829	endif ! of if scavenging
[3203]	1830
[42]	1831	END IF ! of IF (water)
	1832
[2149]	1833	c 9a bis. CO2 clouds (CL & JA)
[1617]	1834	c ---------------------------------------
	1835	c CO2 ice cloud condensation in the atmosphere
	1836	c ----------------------------------------
[1818]	1837	c flag needed in callphys.def:
	1838	c co2clouds=.true. is mandatory (default is .false.)
[3203]	1839	c co2useh2o=.true. if you want to allow co2 condensation
	1840	c on water ice particles
[1818]	1841	c meteo_flux=.true. if you want to add a meteoritic
	1842	c supply of CCN
	1843	c CLFvaryingCO2=.true. if you want to have a sub-grid
[3203]	1844	c temperature distribution
[1818]	1845	c spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti
[3203]	1846	c nuiceco2_sed=0.2 variance of the size distribution for the
[1818]	1847	c sedimentation
[3203]	1848	c nuiceco2_ref=0.2 variance of the size distribution for the
[1818]	1849	c nucleation
	1850	c imicroco2=50 micro-timestep is 1/50 of physical timestep
[2362]	1851	zdqssed_co2(:) = 0.
[2562]	1852	zdqssed_ccn(:,:) = 0.
[1922]	1853
[2362]	1854	IF (co2clouds) THEN
	1855	call co2cloud(ngrid,nlayer,ptimestep,
[1617]	1856	& zplev,zplay,pdpsrf,zzlay,pt,pdt,
	1857	& pq,pdq,zdqcloudco2,zdtcloudco2,
	1858	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
[2447]	1859	& rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2,
[2826]	1860	& zdqssed_ccn,pdu,pu,zcondicea_co2microp)
[1617]	1861
[2562]	1862	DO iq=1, nq
	1863	DO ig=1,ngrid
[2900]	1864	DO islope = 1,nslope
	1865	dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+
	1866	& zdqssed_ccn(ig,iq)cos(pidef_slope_mean(islope)/180.)
	1867	ENDDO !(islope)
[2562]	1868	ENDDO ! (ig)
[2900]	1869	ENDDO ! (iq)q)
[1617]	1870	c Temperature variation due to latent heat release
[2562]	1871	pdt(1:ngrid,1:nlayer) =
[3203]	1872	& pdt(1:ngrid,1:nlayer) +
[2551]	1873	& zdtcloudco2(1:ngrid,1:nlayer)
[1617]	1874
	1875	! increment dust and ccn masses and numbers
	1876	! We need to check that we have Nccn & Ndust > 0
	1877	! This is due to single precision rounding problems
[1816]	1878	! increment dust tracers tendancies
[2362]	1879	pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass)
	1880	& + zdqcloudco2(:,:,igcm_dust_mass)
	1881
	1882	pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number)
	1883	& + zdqcloudco2(:,:,igcm_dust_number)
	1884
	1885	pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2)
	1886	& + zdqcloudco2(:,:,igcm_co2)
	1887
	1888	pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice)
	1889	& + zdqcloudco2(:,:,igcm_co2_ice)
	1890
	1891	pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass)
	1892	& + zdqcloudco2(:,:,igcm_ccnco2_mass)
	1893
	1894	pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number)
	1895	& + zdqcloudco2(:,:,igcm_ccnco2_number)
	1896
[2589]	1897	if (meteo_flux) then
	1898	pdq(:,:,igcm_ccnco2_meteor_mass) =
	1899	& pdq(:,:,igcm_ccnco2_meteor_mass) +
	1900	& zdqcloudco2(:,:,igcm_ccnco2_meteor_mass)
	1901
	1902	pdq(:,:,igcm_ccnco2_meteor_number) =
	1903	& pdq(:,:,igcm_ccnco2_meteor_number)
	1904	& + zdqcloudco2(:,:,igcm_ccnco2_meteor_number)
	1905	end if
[1816]	1906	!Update water ice clouds values as well
	1907	if (co2useh2o) then
[3203]	1908	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
	1909	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
[1816]	1910	& zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
[3203]	1911	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
	1912	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
[1816]	1913	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
[3203]	1914	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
	1915	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
[1816]	1916	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
[2362]	1917
[2562]	1918	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
	1919	& pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
	1920	& zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice)
	1921
	1922	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
	1923	& pdq(:,:,igcm_ccnco2_h2o_mass_ccn) +
	1924	& zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn)
	1925
	1926	pdq(:,:,igcm_ccnco2_h2o_number) =
	1927	& pdq(:,:,igcm_ccnco2_h2o_number) +
	1928	& zdqcloudco2(:,:,igcm_ccnco2_h2o_number)
	1929
[2362]	1930	c Negative values?
	1931	where (pq(:,:,igcm_ccn_mass) +
[1921]	1932	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
[3203]	1933	pdq(:,:,igcm_ccn_mass) =
[1816]	1934	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
[3203]	1935	pdq(:,:,igcm_ccn_number) =
[1816]	1936	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
[1921]	1937	end where
[2362]	1938	c Negative values?
	1939	where (pq(:,:,igcm_ccn_number) +
[1921]	1940	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
[3203]	1941	pdq(:,:,igcm_ccn_mass) =
[1921]	1942	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
[3203]	1943	pdq(:,:,igcm_ccn_number) =
[1921]	1944	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1945	end where
[2562]	1946	where (pq(:,:,igcm_ccnco2_h2o_mass_ice) +
	1947	& pq(:,:,igcm_ccnco2_h2o_mass_ccn) +
	1948	& (pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
	1949	& pdq(:,:,igcm_ccnco2_h2o_mass_ccn)
	1950	& )*ptimestep < 0.)
	1951	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
	1952	& - pq(:,:,igcm_ccnco2_h2o_mass_ice)
	1953	& /ptimestep + 1.e-30
	1954	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
	1955	& - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
	1956	& /ptimestep + 1.e-30
	1957	pdq(:,:,igcm_ccnco2_h2o_number) =
	1958	& - pq(:,:,igcm_ccnco2_h2o_number)
	1959	& /ptimestep + 1.e-30
	1960	end where
	1961
	1962	where (pq(:,:,igcm_ccnco2_h2o_number) +
	1963	& (pdq(:,:,igcm_ccnco2_h2o_number)
	1964	& )*ptimestep < 0.)
	1965	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
	1966	& - pq(:,:,igcm_ccnco2_h2o_mass_ice)
	1967	& /ptimestep + 1.e-30
	1968	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
	1969	& - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
	1970	& /ptimestep + 1.e-30
	1971	pdq(:,:,igcm_ccnco2_h2o_number) =
	1972	& - pq(:,:,igcm_ccnco2_h2o_number)
	1973	& /ptimestep + 1.e-30
	1974	end where
[1921]	1975	endif ! of if (co2useh2o)
[1720]	1976	c Negative values?
[3203]	1977	where (pq(:,:,igcm_ccnco2_mass) +
[1921]	1978	& ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
[3203]	1979	pdq(:,:,igcm_ccnco2_mass) =
[1617]	1980	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
[3203]	1981	pdq(:,:,igcm_ccnco2_number) =
[1617]	1982	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
[1720]	1983	end where
[3203]	1984	where (pq(:,:,igcm_ccnco2_number) +
[1921]	1985	& ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
[3203]	1986	pdq(:,:,igcm_ccnco2_mass) =
[1720]	1987	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
[3203]	1988	pdq(:,:,igcm_ccnco2_number) =
[1720]	1989	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
[1921]	1990	end where
[3203]	1991
[1617]	1992	c Negative values?
[3203]	1993	where (pq(:,:,igcm_dust_mass) +
[1921]	1994	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
[3203]	1995	pdq(:,:,igcm_dust_mass) =
[1922]	1996	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
[3203]	1997	pdq(:,:,igcm_dust_number) =
[1922]	1998	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
[1921]	1999	end where
[3203]	2000	where (pq(:,:,igcm_dust_number) +
[1921]	2001	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
[3203]	2002	pdq(:,:,igcm_dust_mass) =
[1922]	2003	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
[3203]	2004	pdq(:,:,igcm_dust_number) =
[1922]	2005	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
[1921]	2006	end where
[2589]	2007	if (meteo_flux) then
	2008	where (pq(:,:,igcm_ccnco2_meteor_mass) +
	2009	& ptimestep*pdq(:,:,igcm_ccnco2_meteor_mass) < 0.)
	2010	pdq(:,:,igcm_ccnco2_meteor_mass) =
	2011	& - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
	2012	pdq(:,:,igcm_ccnco2_meteor_number) =
	2013	& - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
	2014	end where
	2015	where (pq(:,:,igcm_ccnco2_meteor_number) +
	2016	& ptimestep*pdq(:,:,igcm_ccnco2_meteor_number) < 0.)
	2017	pdq(:,:,igcm_ccnco2_meteor_mass) =
	2018	& - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
	2019	pdq(:,:,igcm_ccnco2_meteor_number) =
	2020	& - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
	2021	end where
	2022	end if
[2362]	2023	END IF ! of IF (co2clouds)
[1617]	2024
[2149]	2025	c 9b. Aerosol particles
[42]	2026	c -------------------
	2027	c ----------
	2028	c Dust devil :
	2029	c ----------
[3203]	2030	IF(callddevil) then
[883]	2031	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
[42]	2032	& zdqdev,zdqsdev)
[2896]	2033
[42]	2034	if (dustbin.ge.1) then
	2035	do iq=1,nq
	2036	DO l=1,nlayer
	2037	DO ig=1,ngrid
	2038	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	2039	ENDDO
	2040	ENDDO
	2041	enddo
	2042	do iq=1,nq
	2043	DO ig=1,ngrid
[2900]	2044	DO islope = 1,nslope
	2045	dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
	2046	& zdqsdev(ig,iq)cos(pidef_slope_mean(islope)/180.)
	2047	ENDDO
[42]	2048	ENDDO
	2049	enddo
	2050	endif ! of if (dustbin.ge.1)
	2051
	2052	END IF ! of IF (callddevil)
	2053
[3203]	2054	c -------------
[42]	2055	c Sedimentation : acts also on water ice
[1974]	2056	c -------------
[3203]	2057	IF (sedimentation) THEN
[42]	2058	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	2059	zdqssed(1:ngrid,1:nq)=0
	2060
[1921]	2061	c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
	2062	c Zdqssed isn't
[2953]	2063
[1974]	2064	call callsedim(ngrid,nlayer,ptimestep,
[2199]	2065	& zplev,zzlev,zzlay,pt,pdt,
	2066	& rdust,rstormdust,rtopdust,
[1974]	2067	& rice,rsedcloud,rhocloud,
[3203]	2068	& pq,pdq,zdqsed,zdqssed,nq,
[411]	2069	& tau,tauscaling)
[2896]	2070
[3157]	2071
[1921]	2072	c Flux at the surface of co2 ice computed in co2cloud microtimestep
[1974]	2073	IF (rdstorm) THEN
	2074	c Storm dust cannot sediment to the surface
[3203]	2075	DO ig=1,ngrid
[1974]	2076	zdqsed(ig,1,igcm_stormdust_mass)=
	2077	& zdqsed(ig,1,igcm_stormdust_mass)+
	2078	& zdqssed(ig,igcm_stormdust_mass) /
	2079	& ((pplev(ig,1)-pplev(ig,2))/g)
	2080	zdqsed(ig,1,igcm_stormdust_number)=
	2081	& zdqsed(ig,1,igcm_stormdust_number)+
	2082	& zdqssed(ig,igcm_stormdust_number) /
[3203]	2083	& ((pplev(ig,1)-pplev(ig,2))/g)
[1974]	2084	zdqssed(ig,igcm_stormdust_mass)=0.
	2085	zdqssed(ig,igcm_stormdust_number)=0.
	2086	ENDDO
	2087	ENDIF !rdstorm
	2088
[42]	2089	DO iq=1, nq
	2090	DO l=1,nlayer
	2091	DO ig=1,ngrid
	2092	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	2093	ENDDO
	2094	ENDDO
	2095	ENDDO
	2096	DO iq=1, nq
	2097	DO ig=1,ngrid
[2900]	2098	DO islope = 1,nslope
	2099	dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
	2100	& zdqssed(ig,iq)cos(pidef_slope_mean(islope)/180.)
	2101	ENDDO
[42]	2102	ENDDO
	2103	ENDDO
[1974]	2104
[42]	2105	END IF ! of IF (sedimentation)
[1974]	2106
[1236]	2107	c Add lifted dust to tendancies after sedimentation in the LES (AC)
	2108	IF (turb_resolved) THEN
	2109	DO iq=1, nq
	2110	DO l=1,nlayer
	2111	DO ig=1,ngrid
	2112	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	2113	ENDDO
	2114	ENDDO
	2115	ENDDO
	2116	DO iq=1, nq
	2117	DO ig=1,ngrid
[2900]	2118	dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
[1236]	2119	ENDDO
	2120	ENDDO
	2121	ENDIF
[556]	2122	c
[2149]	2123	c 9c. Chemical species
[556]	2124	c ------------------
[42]	2125
[556]	2126	#ifndef MESOSCALE
	2127	c --------------
	2128	c photochemistry :
	2129	c --------------
[2158]	2130	IF (photochem) then
[411]	2131
[2164]	2132	if (modulo(icount-1,ichemistry).eq.0) then
	2133	! compute chemistry every ichemistry physics step
	2134
[556]	2135	! dust and ice surface area
[1246]	2136	call surfacearea(ngrid, nlayer, naerkind,
[3203]	2137	$ ptimestep, zplay, zzlay,
	2138	$ pt, pq, pdq, nq,
	2139	$ rdust, rice, tau, tauscaling,
[556]	2140	$ surfdust, surfice)
	2141	! call photochemistry
[2900]	2142	DO ig = 1,ngrid
	2143	qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig))
	2144	ENDDO
[1047]	2145	call calchim(ngrid,nlayer,nq,
[1036]	2146	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
[556]	2147	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
[2900]	2148	$ zdqcloud,zdqscloud,tau(:,1),
	2149	$ qsurf_tmp(:,igcm_co2),
[556]	2150	$ pu,pdu,pv,pdv,surfdust,surfice)
[2896]	2151
[2164]	2152	endif ! of if (modulo(icount-1,ichemistry).eq.0)
	2153
[556]	2154	! increment values of tracers:
	2155	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	2156	! tracers is zero anyways
	2157	DO l=1,nlayer
	2158	DO ig=1,ngrid
	2159	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	2160	ENDDO
	2161	ENDDO
	2162	ENDDO ! of DO iq=1,nq
[3203]	2163
[556]	2164	! add condensation tendency for H2O2
	2165	if (igcm_h2o2.ne.0) then
	2166	DO l=1,nlayer
	2167	DO ig=1,ngrid
	2168	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	2169	& +zdqcloud(ig,l,igcm_h2o2)
	2170	ENDDO
	2171	ENDDO
	2172	endif
	2173
	2174	! increment surface values of tracers:
	2175	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	2176	! tracers is zero anyways
	2177	DO ig=1,ngrid
[2900]	2178	DO islope = 1,nslope
	2179	dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) +
	2180	& zdqschim(ig,iq)cos(pidef_slope_mean(islope)/180.)
	2181	ENDDO
[556]	2182	ENDDO
	2183	ENDDO ! of DO iq=1,nq
	2184
	2185	! add condensation tendency for H2O2
	2186	if (igcm_h2o2.ne.0) then
	2187	DO ig=1,ngrid
[2900]	2188	DO islope = 1,nslope
	2189	dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+
	2190	& zdqscloud(ig,igcm_h2o2)cos(pidef_slope_mean(islope)/180.)
	2191	ENDDO
[556]	2192	ENDDO
	2193	endif
	2194
[2158]	2195	END IF ! of IF (photochem)
[556]	2196	#endif
	2197
[42]	2198
[226]	2199	#ifndef MESOSCALE
	2200	c-----------------------------------------------------------------------
[2149]	2201	c 10. THERMOSPHERE CALCULATION
[226]	2202	c-----------------------------------------------------------------------
[42]	2203
[226]	2204	if (callthermos) then
[1047]	2205	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
[2908]	2206	$ mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay,
[226]	2207	& pt,pq,pu,pv,pdt,pdq,
[2467]	2208	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
	2209	$ PhiEscH,PhiEscH2,PhiEscD)
[226]	2210
	2211	DO l=1,nlayer
	2212	DO ig=1,ngrid
	2213	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
[2362]	2214	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
[226]	2215	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	2216	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	2217	DO iq=1, nq
	2218	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	2219	ENDDO
	2220	ENDDO
	2221	ENDDO
	2222
	2223	endif ! of if (callthermos)
	2224	#endif
[2963]	2225
[42]	2226	c-----------------------------------------------------------------------
[2149]	2227	c 11. Carbon dioxide condensation-sublimation:
[1114]	2228	c (should be the last atmospherical physical process to be computed)
	2229	c -------------------------------------------
	2230	IF (tituscap) THEN
	2231	!!! get the actual co2 seasonal cap from Titus observations
[2362]	2232	CALL geticecover(ngrid, 180.*zls/pi,
[2826]	2233	. 180.longitude/pi, 180.latitude/pi,
[2900]	2234	. qsurf_tmp(:,igcm_co2) )
	2235	qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
[1114]	2236	ENDIF
	2237
[3203]	2238
[1114]	2239	IF (callcond) THEN
[2362]	2240	zdtc(:,:) = 0.
[2900]	2241	zdtsurfc(:,:) = 0.
[2362]	2242	zduc(:,:) = 0.
	2243	zdvc(:,:) = 0.
	2244	zdqc(:,:,:) = 0.
[2900]	2245	zdqsc(:,:,:) = 0.
[2953]	2246	CALL co2condens(ngrid,nlayer,nq,nslope,ptimestep,
[1114]	2247	$ capcal,zplay,zplev,tsurf,pt,
	2248	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
[3130]	2249	$ qsurf(:,igcm_co2,:),perennial_co2ice,
[2999]	2250	$ albedo,emis,rdust,
[1114]	2251	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
[2685]	2252	$ fluxsurf_dn_sw,zls,
[1996]	2253	$ zdqssed_co2,zcondicea_co2microp,
[2551]	2254	& zdqsc)
[2896]	2255
[3207]	2256	if (ngrid == 1) then ! For the 1D model
	2257	! CO2cond_ps is a coefficient to control the surface pressure change
	2258	pdpsrf = CO2cond_ps*pdpsrf
	2259	zduc = CO2cond_ps*zduc
	2260	zdvc = CO2cond_ps*zdvc
	2261	zdqc = CO2cond_ps*zdqc
	2262	endif
	2263
[2566]	2264	DO iq=1, nq
[1114]	2265	DO ig=1,ngrid
[2900]	2266	dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:)
[2362]	2267	ENDDO ! (ig)
[2566]	2268	ENDDO ! (iq)
	2269	DO l=1,nlayer
[2362]	2270	DO ig=1,ngrid
[1114]	2271	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	2272	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	2273	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	2274	ENDDO
[2566]	2275	ENDDO
	2276	DO ig=1,ngrid
[2900]	2277	zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:)
[2566]	2278	ENDDO
[1114]	2279
	2280	DO iq=1, nq
	2281	DO l=1,nlayer
	2282	DO ig=1,ngrid
[2362]	2283	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
[1114]	2284	ENDDO
	2285	ENDDO
	2286	ENDDO
	2287
	2288	#ifndef MESOSCALE
	2289	! update surface pressure
	2290	DO ig=1,ngrid
	2291	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
	2292	ENDDO
	2293	! update pressure levels
	2294	DO l=1,nlayer
	2295	DO ig=1,ngrid
	2296	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
	2297	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
	2298	ENDDO
	2299	ENDDO
	2300	zplev(:,nlayer+1) = 0.
[3157]	2301
[3203]	2302	! Calculation of zzlay and zzlay with udpated pressure and temperature
	2303	DO ig=1,ngrid
	2304	zzlay(ig,1)=-(log(zplay(ig,1)/ps(ig)))rnew(ig,1)
	2305	& (pt(ig,1)+pdt(ig,1)*ptimestep) /g
	2306
	2307	DO l=2,nlayer
	2308
[3157]	2309	! compute "mean" temperature of the layer
	2310	if((pt(ig,l)+pdt(ig,l)*ptimestep) .eq.
	2311	& (pt(ig,l-1)+pdt(ig,l-1)*ptimestep)) then
	2312	tlaymean= pt(ig,l)+pdt(ig,l)*ptimestep
	2313	else
[3203]	2314	tlaymean=((pt(ig,l)+pdt(ig,l)*ptimestep)-
	2315	& (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))/
	2316	& log((pt(ig,l)+pdt(ig,l)*ptimestep)/
	2317	& (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))
[3157]	2318	endif
[3203]	2319
[3157]	2320	! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
	2321	gz(ig,l)=g(rad2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))*2
	2322
[3203]	2323
	2324	zzlay(ig,l)=zzlay(ig,l-1)-
[3157]	2325	& (log(zplay(ig,l)/zplay(ig,l-1))rnew(ig,l)tlaymean/gz(ig,l))
[3203]	2326
	2327
[3157]	2328	! update layers altitude
	2329	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
	2330	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
	2331	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
[1114]	2332	ENDDO
[3203]	2333	ENDDO
[1114]	2334	#endif
[1996]	2335	ENDIF ! of IF (callcond)
	2336
[1114]	2337	c-----------------------------------------------------------------------
[2184]	2338	c Updating tracer budget on surface
[3203]	2339	c-----------------------------------------------------------------------
[2184]	2340	DO iq=1, nq
	2341	DO ig=1,ngrid
[2900]	2342	DO islope = 1,nslope
	2343	qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
	2344	& ptimestep*dqsurf(ig,iq,islope)
[3203]	2345	ENDDO
[2184]	2346	ENDDO ! (ig)
	2347	ENDDO ! (iq)
	2348	c-----------------------------------------------------------------------
[2149]	2349	c 12. Surface and sub-surface soil temperature
[42]	2350	c-----------------------------------------------------------------------
	2351	c
	2352	c
[2149]	2353	c 12.1 Increment Surface temperature:
[42]	2354	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2355
	2356	DO ig=1,ngrid
[2900]	2357	DO islope = 1,nslope
	2358	tsurf(ig,islope)=tsurf(ig,islope)+
[3203]	2359	& ptimestep*zdtsurf(ig,islope)
[2900]	2360	ENDDO
[42]	2361	ENDDO
	2362
	2363	c Prescribe a cold trap at south pole (except at high obliquity !!)
	2364	c Temperature at the surface is set there to be the temperature
	2365	c corresponding to equilibrium temperature between phases of CO2
	2366
[411]	2367
[2823]	2368	IF (water) THEN
[3203]	2369	!#ifndef MESOSCALE
[2009]	2370	! if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
[226]	2371	! NB: Updated surface pressure, at grid point 'ngrid', is
[883]	2372	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
[1114]	2373	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	2374	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	2375	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
	2376	! endif
	2377	!#endif
[42]	2378	c -------------------------------------------------------------
[283]	2379	c Change of surface albedo in case of ground frost
[42]	2380	c everywhere except on the north permanent cap and in regions
[3203]	2381	c covered by dry ice.
[2009]	2382	c ALWAYS PLACE these lines after co2condens !!!
[42]	2383	c -------------------------------------------------------------
[3586]	2384	do ig = 1,ngrid
[2900]	2385	do islope = 1,nslope
[3586]	2386	if (abs(qsurf(ig,igcm_co2,islope)) < 1.e-10) then ! No CO2 frost
	2387
	2388	if (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)
	2389	& then ! There is H2O frost
	2390	if (cst_cap_albedo .and. watercaptag(ig) .and.
	2391	& abs(perennial_co2ice(ig,islope)) < 1.e-10) then ! Water cap remains unchanged by water frost deposition and no CO2 perennial ice
	2392	albedo(ig,:,islope) = albedo_h2o_cap
	2393	emis(ig,islope) = 1.
	2394	else
	2395	albedo(ig,:,islope) = albedo_h2o_frost
	2396	emis(ig,islope) = 1.
	2397	endif
	2398	else ! No H2O frost
	2399	if (abs(perennial_co2ice(ig,islope)) < 1.e-10 .and.
	2400	& watercaptag(ig)) then ! No CO2 perennial ice but there is water cap
	2401	albedo(ig,:,islope) = albedo_h2o_cap
	2402	emis(ig,islope) = 1.
	2403	endif
	2404	endif
	2405
[42]	2406	endif
[2900]	2407	enddo ! islope
[42]	2408	enddo ! of do ig=1,ngrid
[2823]	2409	ENDIF ! of IF (water)
[42]	2410
[528]	2411	c
[2149]	2412	c 12.2 Compute soil temperatures and subsurface heat flux:
[42]	2413	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2414	IF (callsoil) THEN
[833]	2415	c Thermal inertia feedback
[3230]	2416	IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
	2417
	2418	CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
	2419	s qsurf(:,igcm_h2o_ice,:),pore_icefraction,
	2420	s inertiesoil_tifeedback(:,:,:))
	2421
[2942]	2422	CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback,
[833]	2423	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	2424	ELSE
[2942]	2425	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
[833]	2426	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	2427	ENDIF
[42]	2428	ENDIF
	2429
[1974]	2430	c To avoid negative values
	2431	IF (rdstorm) THEN
[3203]	2432	where (pq(:,:,igcm_stormdust_mass) +
[1974]	2433	& ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
[3203]	2434	pdq(:,:,igcm_stormdust_mass) =
[1974]	2435	& - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
[3203]	2436	pdq(:,:,igcm_stormdust_number) =
[1974]	2437	& - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
	2438	end where
[3203]	2439	where (pq(:,:,igcm_stormdust_number) +
[1974]	2440	& ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
[3203]	2441	pdq(:,:,igcm_stormdust_mass) =
[1974]	2442	& - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
[3203]	2443	pdq(:,:,igcm_stormdust_number) =
[1974]	2444	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	2445	end where
	2446
[3203]	2447	where (pq(:,:,igcm_dust_mass) +
[1974]	2448	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
[3203]	2449	pdq(:,:,igcm_dust_mass) =
[1974]	2450	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
[3203]	2451	pdq(:,:,igcm_dust_number) =
[1974]	2452	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	2453	end where
[3203]	2454	where (pq(:,:,igcm_dust_number) +
[1974]	2455	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
[3203]	2456	pdq(:,:,igcm_dust_mass) =
[1974]	2457	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
[3203]	2458	pdq(:,:,igcm_dust_number) =
[1974]	2459	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	2460	end where
[3203]	2461	ENDIF !(rdstorm)
[1974]	2462
[42]	2463	c-----------------------------------------------------------------------
[2260]	2464	c J. Naar : Surface and sub-surface water ice
	2465	c-----------------------------------------------------------------------
	2466	c
	2467	c
	2468	c Increment Watercap (surface h2o reservoirs):
	2469	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2470
	2471	DO ig=1,ngrid
[2900]	2472	DO islope = 1,nslope
	2473	watercap(ig,islope)=watercap(ig,islope)+
	2474	s ptimestep*dwatercap(ig,islope)
	2475	ENDDO
[2260]	2476	ENDDO
	2477
[2561]	2478	IF (refill_watercap) THEN
	2479
[3159]	2480	DO ig = 1,ngrid
[2900]	2481	DO islope = 1,nslope
[3159]	2482	if (watercaptag(ig) .and. (qsurf(ig,igcm_h2o_ice,islope)
	2483	& > frost_metam_threshold)) then
[2561]	2484
[3159]	2485	watercap(ig,islope) = watercap(ig,islope)
	2486	& + qsurf(ig,igcm_h2o_ice,islope)
	2487	& - frost_metam_threshold
[2900]	2488	qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold
[3159]	2489	endif ! watercaptag
[2900]	2490	ENDDO
[2561]	2491	ENDDO
	2492
[3159]	2493	ENDIF ! refill_watercap
[2561]	2494
[2260]	2495	c-----------------------------------------------------------------------
[2149]	2496	c 13. Write output files
[42]	2497	c ----------------------
[2900]	2498	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
	2499	& albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
	2500
[42]	2501	c -------------------------------
	2502	c Dynamical fields incrementation
	2503	c -------------------------------
	2504	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	2505	! temperature, zonal and meridional wind
	2506	DO l=1,nlayer
	2507	DO ig=1,ngrid
	2508	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	2509	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	2510	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	2511	ENDDO
	2512	ENDDO
	2513
	2514	! tracers
	2515	DO iq=1, nq
	2516	DO l=1,nlayer
	2517	DO ig=1,ngrid
	2518	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	2519	ENDDO
	2520	ENDDO
	2521	ENDDO
	2522
[2362]	2523	! Density
[42]	2524	DO l=1,nlayer
	2525	DO ig=1,ngrid
	2526	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	2527	ENDDO
	2528	ENDDO
	2529
[269]	2530	! Potential Temperature
	2531
[1047]	2532	DO ig=1,ngrid
	2533	DO l=1,nlayer
[528]	2534	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
[269]	2535	ENDDO
	2536	ENDDO
	2537
[226]	2538	c Compute surface stress : (NB: z0 is a common in surfdat.h)
[42]	2539	c DO ig=1,ngrid
[226]	2540	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
[42]	2541	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	2542	c ENDDO
	2543
	2544	c Sum of fluxes in solar spectral bands (for output only)
[2685]	2545	fluxtop_dn_sw_tot(1:ngrid)=fluxtop_dn_sw(1:ngrid,1) +
	2546	& fluxtop_dn_sw(1:ngrid,2)
	2547	fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) +
	2548	& fluxtop_up_sw(1:ngrid,2)
[2900]	2549	fluxsurf_dn_sw_tot(1:ngrid,1:nslope)=
	2550	& fluxsurf_dn_sw(1:ngrid,1,1:nslope) +
	2551	& fluxsurf_dn_sw(1:ngrid,2,1:nslope)
[2685]	2552	fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) +
	2553	& fluxsurf_up_sw(1:ngrid,2)
	2554
[42]	2555	c ***** TEST ****************************************************
[835]	2556	ztim1 = 999
	2557	DO l=1,nlayer
	2558	DO ig=1,ngrid
	2559	if (pt(ig,l).lt.ztim1) then
	2560	ztim1 = pt(ig,l)
	2561	igmin = ig
[3203]	2562	lmin = l
[835]	2563	end if
	2564	ENDDO
	2565	ENDDO
	2566	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
	2567	write(,) 'PHYSIQ: stability WARNING :'
	2568	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	2569	& 'ig l =', igmin, lmin
	2570	end if
[42]	2571
[566]	2572	c ----------------------------------------------------------
	2573	c ----------------------------------------------------------
	2574	c INTERPOLATIONS IN THE SURFACE-LAYER
	2575	c ----------------------------------------------------------
	2576	c ----------------------------------------------------------
	2577
[657]	2578	n_out=0 ! number of elements in the z_out array.
[636]	2579	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
	2580	! to 5
	2581	IF (n_out .ne. 0) THEN
[566]	2582
[636]	2583	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
	2584	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
	2585	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
[566]	2586
[636]	2587	z_out(:)=[3.,2.,1.,0.5,0.1]
	2588	u_out(:,:)=0.
	2589	T_out(:,:)=0.
[566]	2590
[636]	2591	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
[3219]	2592	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,q2,tsurf(:,iflat),
[3325]	2593	& zh,zq(:,1,igcm_h2o_vap),qsurf(:,igcm_h2o_ice,iflat),mmean(:,1),
	2594	& z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv)
[1242]	2595	! pourquoi ustar recalcule ici? fait dans vdifc.
[566]	2596
	2597	#ifndef MESOSCALE
	2598	DO n=1,n_out
[636]	2599	write(zstring, '(F8.6)') z_out(n)
[2932]	2600	call write_output('T_out_'//trim(zstring),
	2601	& 'potential temperature at z_out','K',T_out(:,n))
	2602	call write_output('u_out_'//trim(zstring),
	2603	& 'horizontal velocity norm at z_out','m/s',u_out(:,n))
[636]	2604	ENDDO
[2932]	2605	call write_output('u_star',
	2606	& 'friction velocity','m/s',ustar)
	2607	call write_output('teta_star',
	2608	& 'friction potential temperature','K',tstar)
	2609	call write_output('vvv',
	2610	& 'Vertical velocity variance at zout','m',vvv)
	2611	call write_output('vhf',
	2612	& 'Vertical heat flux at zout','m',vhf)
[636]	2613	#else
	2614	T_out1(:)=T_out(:,1)
	2615	u_out1(:)=u_out(:,1)
[566]	2616	#endif
	2617
	2618	ENDIF
	2619
	2620	c ----------------------------------------------------------
	2621	c ----------------------------------------------------------
	2622	c END OF SURFACE LAYER INTERPOLATIONS
	2623	c ----------------------------------------------------------
	2624	c ----------------------------------------------------------
	2625
[226]	2626	#ifndef MESOSCALE
[42]	2627	c -------------------------------------------------------------------
	2628	c Writing NetCDF file "RESTARTFI" at the end of the run
	2629	c -------------------------------------------------------------------
	2630	c Note: 'restartfi' is stored just before dynamics are stored
	2631	c in 'restart'. Between now and the writting of 'restart',
	2632	c there will have been the itau=itau+1 instruction and
	2633	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	2634	c thus we store for time=time+dtvr
	2635
[2646]	2636	! default: not writing a restart file at this time step
[3203]	2637	write_restart=.false.
[2646]	2638	IF (ecritstart.GT.0) THEN
	2639	! For when we store multiple time steps in the restart file
	2640	IF (MODULO(icount*iphysiq,ecritstart).EQ.0) THEN
	2641	write_restart=.true.
	2642	ENDIF
	2643	ENDIF
	2644	IF (lastcall) THEN
	2645	! Always write a restart at the end of the simulation
	2646	write_restart=.true.
	2647	ENDIF
[3203]	2648
[2646]	2649	IF (write_restart) THEN
	2650	IF (grid_type==unstructured) THEN !IF DYNAMICO
[2507]	2651
[2501]	2652	! When running Dynamico, no need to add a dynamics time step to ztime_fin
[3203]	2653	IF (ptime.LE. 1.E-10) THEN
[2501]	2654	! Residual ptime occurs with Dynamico
	2655	ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec)
	2656	. - day_ini - time_phys
	2657	ELSE
	2658	ztime_fin = pday + ptime !+ ptimestep/(float(iphysiq)*daysec)
	2659	. - day_ini - time_phys
	2660	ENDIF
[2514]	2661	if (ecritstart==0) then
	2662	ztime_fin = ztime_fin-(day_end-day_ini)
	2663	endif
[2507]	2664
[2646]	2665	ELSE ! IF LMDZ
[2507]	2666
[3094]	2667	if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
[3203]	2668	ztime_fin = pday - day_ini + ptime
[3094]	2669	& + ptimestep/(float(iphysiq)*daysec)
	2670	else !IF ONE RESTART final time in top of day_end
[3203]	2671	ztime_fin = pday - day_ini-(day_end-day_ini)
[3094]	2672	& + ptime + ptimestep/(float(iphysiq)*daysec)
	2673	endif
[2507]	2674
[2646]	2675	ENDIF ! of IF (grid_type==unstructured)
[3203]	2676	write(*,'(A,I7,A,F12.5)')
[2646]	2677	. 'PHYSIQ: writing in restartfi ; icount=',
[999]	2678	. icount,' date=',ztime_fin
[3203]	2679
[3113]	2680	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,
[999]	2681	. ptimestep,ztime_fin,
[2942]	2682	. tsurf,tsoil,inertiesoil,albedo,
[3113]	2683	. emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,
[3130]	2684	. wstar,watercap,perennial_co2ice)
[2924]	2685	ENDIF ! of IF (write_restart)
	2686
[226]	2687	#endif
[42]	2688
[2931]	2689	c IF (ngrid.NE.1) then
[2924]	2690
[42]	2691	c -------------------------------------------------------------------
	2692	c Calculation of diagnostic variables written in both stats and
	2693	c diagfi files
	2694	c -------------------------------------------------------------------
[2678]	2695	do ig=1,ngrid
[3203]	2696	if(mu0(ig).le.0.01) then
	2697	fluxsurf_dir_dn_sw(ig) = 0.
	2698	else
[3147]	2699	if (water) then
	2700	! both water and dust contribute
[3148]	2701	fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_solmu0(ig)
[2678]	2702	& exp(-(tau(ig,iaer_dust_doubleq)+
	2703	& tau(ig,iaer_h2o_ice))/mu0(ig))
[3147]	2704	else
	2705	! only dust contributes
[3148]	2706	fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_solmu0(ig)
[3147]	2707	& exp(-(tau(ig,iaer_dust_doubleq))/mu0(ig))
	2708	endif ! of if (water)
	2709	endif ! of if(mu0(ig).le.0.01)
[2678]	2710	enddo
	2711
[2282]	2712	! Density-scaled opacities
[2362]	2713	do ig=1,ngrid
	2714	dsodust(ig,:) =
[2282]	2715	& dsodust(ig,:)*tauscaling(ig)
[2362]	2716	dsords(ig,:) =
[2282]	2717	& dsords(ig,:)*tauscaling(ig)
	2718	dsotop(ig,:) =
	2719	& dsotop(ig,:)*tauscaling(ig)
	2720	enddo
[2362]	2721
[719]	2722	if(doubleq) then
[3203]	2723	do ig=1,ngrid
	2724	IF (sedimentation) THEN
	2725	dqdustsurf(ig) =
[756]	2726	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
[3203]	2727	dndustsurf(ig) =
[756]	2728	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
[2931]	2729	ENDIF
[756]	2730	ndust(ig,:) =
[1720]	2731	& zq(ig,:,igcm_dust_number)*tauscaling(ig)
[756]	2732	qdust(ig,:) =
[1720]	2733	& zq(ig,:,igcm_dust_mass)*tauscaling(ig)
[719]	2734	enddo
	2735	if (scavenging) then
[3203]	2736	do ig=1,ngrid
	2737	IF (sedimentation) THEN
	2738	dqdustsurf(ig) = dqdustsurf(ig) +
[756]	2739	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
[3203]	2740	dndustsurf(ig) = dndustsurf(ig) +
[756]	2741	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
[2931]	2742	ENDIF
[756]	2743	nccn(ig,:) =
[1720]	2744	& zq(ig,:,igcm_ccn_number)*tauscaling(ig)
[756]	2745	qccn(ig,:) =
[1720]	2746	& zq(ig,:,igcm_ccn_mass)*tauscaling(ig)
[719]	2747	enddo
	2748	endif
[1921]	2749	endif ! of (doubleq)
[1974]	2750
	2751	if (rdstorm) then ! diagnostics of stormdust tendancies for 1D and 3D
	2752	mstormdtot(:)=0
	2753	mdusttot(:)=0
	2754	qdusttotal(:,:)=0
[3203]	2755	do ig=1,ngrid
	2756	rdsdqdustsurf(ig) =
[1974]	2757	& zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
[3203]	2758	rdsdndustsurf(ig) =
[1974]	2759	& zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
	2760	rdsndust(ig,:) =
	2761	& pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
	2762	rdsqdust(ig,:) =
	2763	& pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
	2764	do l=1,nlayer
[3203]	2765	mstormdtot(ig) = mstormdtot(ig) +
	2766	& zq(ig,l,igcm_stormdust_mass) *
[1974]	2767	& (zplev(ig,l) - zplev(ig,l+1)) / g
[3203]	2768	mdusttot(ig) = mdusttot(ig) +
	2769	& zq(ig,l,igcm_dust_mass) *
[1974]	2770	& (zplev(ig,l) - zplev(ig,l+1)) / g
	2771	qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
	2772	enddo
	2773	enddo
	2774	endif !(rdstorm)
[3203]	2775
[42]	2776	if (water) then
[286]	2777	mtot(:)=0
	2778	icetot(:)=0
	2779	rave(:)=0
	2780	tauTES(:)=0
[2312]	2781
	2782	IF (hdo) then
	2783	mtotD(:)=0
	2784	icetotD(:)=0
	2785	ENDIF !hdo
	2786
[2362]	2787	do ig=1,ngrid
[1047]	2788	do l=1,nlayer
[3203]	2789	mtot(ig) = mtot(ig) +
	2790	& zq(ig,l,igcm_h2o_vap) *
[883]	2791	& (zplev(ig,l) - zplev(ig,l+1)) / g
[3203]	2792	icetot(ig) = icetot(ig) +
	2793	& zq(ig,l,igcm_h2o_ice) *
[883]	2794	& (zplev(ig,l) - zplev(ig,l+1)) / g
[2312]	2795	IF (hdo) then
[2362]	2796	mtotD(ig) = mtotD(ig) +
	2797	& zq(ig,l,igcm_hdo_vap) *
[2312]	2798	& (zplev(ig,l) - zplev(ig,l+1)) / g
[2362]	2799	icetotD(ig) = icetotD(ig) +
	2800	& zq(ig,l,igcm_hdo_ice) *
[2312]	2801	& (zplev(ig,l) - zplev(ig,l+1)) / g
	2802	ENDIF !hdo
	2803
[42]	2804	c Computing abs optical depth at 825 cm-1 in each
	2805	c layer to simulate NEW TES retrieval
	2806	Qabsice = min(
	2807	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	2808	& )
[3203]	2809	opTES(ig,l)= 0.75 * Qabsice *
[42]	2810	& zq(ig,l,igcm_h2o_ice) *
[883]	2811	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	2812	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
[3203]	2813	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
[42]	2814	enddo
[669]	2815	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
	2816	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
[42]	2817	enddo
[1047]	2818	call watersat(ngrid*nlayer,zt,zplay,zqsat)
[669]	2819	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
[42]	2820
[669]	2821	if (scavenging) then
	2822	Nccntot(:)= 0
	2823	Mccntot(:)= 0
	2824	rave(:)=0
[3203]	2825	do ig=1,ngrid
[1047]	2826	do l=1,nlayer
[3203]	2827	Nccntot(ig) = Nccntot(ig) +
[669]	2828	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
[883]	2829	& *(zplev(ig,l) - zplev(ig,l+1)) / g
[3203]	2830	Mccntot(ig) = Mccntot(ig) +
[669]	2831	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
[883]	2832	& *(zplev(ig,l) - zplev(ig,l+1)) / g
[3203]	2833	cccc Column integrated effective ice radius
	2834	cccc is weighted by total ice surface area (BETTER than total ice mass)
	2835	rave(ig) = rave(ig) +
[669]	2836	& tauscaling(ig) *
	2837	& zq(ig,l,igcm_ccn_number) *
[3203]	2838	& (zplev(ig,l) - zplev(ig,l+1)) / g *
[669]	2839	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
	2840	enddo
	2841	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
	2842	& /max(pi*rave(ig),1.e-30) ! surface weight
	2843	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	2844	enddo
[833]	2845	else ! of if (scavenging)
	2846	rave(:)=0
[3203]	2847	do ig=1,ngrid
[1047]	2848	do l=1,nlayer
[3203]	2849	rave(ig) = rave(ig) +
[833]	2850	& zq(ig,l,igcm_h2o_ice) *
[3203]	2851	& (zplev(ig,l) - zplev(ig,l+1)) / g *
[833]	2852	& rice(ig,l) * (1.+nuice_ref)
	2853	enddo
[3203]	2854	rave(ig) = max(rave(ig) /
[833]	2855	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
	2856	enddo
[669]	2857	endif ! of if (scavenging)
	2858
[1711]	2859	!Alternative A. Pottier weighting
	2860	rave2(:) = 0.
	2861	totrave2(:) = 0.
	2862	do ig=1,ngrid
	2863	do l=1,nlayer
	2864	rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
	2865	totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
	2866	end do
	2867	rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
	2868	end do
	2869
[42]	2870	endif ! of if (water)
[2551]	2871
	2872	if (co2clouds) then
	2873	mtotco2(1:ngrid) = 0.
	2874	icetotco2(1:ngrid) = 0.
	2875	vaptotco2(1:ngrid) = 0.
	2876	do ig=1,ngrid
	2877	do l=1,nlayer
[3203]	2878	vaptotco2(ig) = vaptotco2(ig) +
	2879	& zq(ig,l,igcm_co2) *
[2551]	2880	& (zplev(ig,l) - zplev(ig,l+1)) / g
[3203]	2881	icetotco2(ig) = icetot(ig) +
	2882	& zq(ig,l,igcm_co2_ice) *
[2551]	2883	& (zplev(ig,l) - zplev(ig,l+1)) / g
	2884	end do
	2885	mtotco2(ig) = icetotco2(ig) + vaptotco2(ig)
	2886	end do
	2887	end if
[2823]	2888
[1212]	2889	#ifndef MESOSCALE
[42]	2890	c -----------------------------------------------------------------
	2891	c WSTATS: Saving statistics
	2892	c -----------------------------------------------------------------
[2563]	2893	c ("stats" stores and accumulates key variables in file "stats.nc"
[42]	2894	c which can later be used to make the statistic files of the run:
[2563]	2895	c if flag "callstats" from callphys.def is .true.)
[3203]	2896
[695]	2897	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
[2900]	2898	call wstats(ngrid,"tsurf","Surface temperature","K",2
	2899	& ,tsurf(:,iflat))
[695]	2900	call wstats(ngrid,"co2ice","CO2 ice cover",
[2900]	2901	& "kg.m-2",2,qsurf(:,igcm_co2,iflat))
[2260]	2902	call wstats(ngrid,"watercap","H2O ice cover",
[2900]	2903	& "kg.m-2",2,watercap(:,iflat))
[2415]	2904	call wstats(ngrid,"tau_pref_scenario",
	2905	& "prescribed visible dod at 610 Pa","NU",
	2906	& 2,tau_pref_scenario)
	2907	call wstats(ngrid,"tau_pref_gcm",
	2908	& "visible dod at 610 Pa in the GCM","NU",
	2909	& 2,tau_pref_gcm)
[695]	2910	call wstats(ngrid,"fluxsurf_lw",
[575]	2911	& "Thermal IR radiative flux to surface","W.m-2",2,
[2900]	2912	& fluxsurf_lw(:,iflat))
[2685]	2913	call wstats(ngrid,"fluxsurf_dn_sw",
	2914	& "Incoming Solar radiative flux to surface","W.m-2",2,
[2900]	2915	& fluxsurf_dn_sw_tot(:,iflat))
[2685]	2916	call wstats(ngrid,"fluxsurf_up_sw",
	2917	& "Reflected Solar radiative flux from surface","W.m-2",2,
	2918	& fluxsurf_up_sw_tot)
[695]	2919	call wstats(ngrid,"fluxtop_lw",
[575]	2920	& "Thermal IR radiative flux to space","W.m-2",2,
	2921	& fluxtop_lw)
[2685]	2922	call wstats(ngrid,"fluxtop_dn_sw",
	2923	& "Incoming Solar radiative flux from space","W.m-2",2,
	2924	& fluxtop_dn_sw_tot)
	2925	call wstats(ngrid,"fluxtop_up_sw",
	2926	& "Outgoing Solar radiative flux to space","W.m-2",2,
	2927	& fluxtop_up_sw_tot)
[695]	2928	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	2929	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	2930	call wstats(ngrid,"v","Meridional (North-South) wind",
[575]	2931	& "m.s-1",3,zv)
[705]	2932	call wstats(ngrid,"w","Vertical (down-up) wind",
	2933	& "m.s-1",3,pw)
[695]	2934	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
[883]	2935	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
[758]	2936	call wstats(ngrid,"q2",
	2937	& "Boundary layer eddy kinetic energy",
	2938	& "m2.s-2",3,q2)
	2939	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
[2900]	2940	& emis(:,iflat))
[2685]	2941	call wstats(ngrid,"fluxsurf_dir_dn_sw",
	2942	& "Direct incoming SW flux at surface",
[3203]	2943	& "W.m-2",2,fluxsurf_dir_dn_sw)
[42]	2944
[758]	2945	if (calltherm) then
	2946	call wstats(ngrid,"zmax_th","Height of thermals",
	2947	& "m",2,zmax_th)
	2948	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
	2949	& "K.m/s",2,hfmax_th)
	2950	call wstats(ngrid,"wstar",
	2951	& "Max vertical velocity in thermals",
	2952	& "m/s",2,wstar)
	2953	endif
	2954
[226]	2955	if (water) then
[1047]	2956	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
	2957	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
[719]	2958	call wstats(ngrid,"vmr_h2ovap",
[520]	2959	& "H2O vapor volume mixing ratio","mol/mol",
	2960	& 3,vmr)
[1047]	2961	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	2962	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
[520]	2963	call wstats(ngrid,"vmr_h2oice",
	2964	& "H2O ice volume mixing ratio","mol/mol",
	2965	& 3,vmr)
[1278]	2966	! also store vmr_ice*rice for better diagnostics of rice
	2967	vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
[1502]	2968	& rice(1:ngrid,1:nlayer)
[1278]	2969	call wstats(ngrid,"vmr_h2oice_rice",
	2970	& "H2O ice mixing ratio times ice particule size",
	2971	& "(mol/mol)*m",
	2972	& 3,vmr)
[1047]	2973	vmr=zqsat(1:ngrid,1:nlayer)
	2974	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
[669]	2975	call wstats(ngrid,"vmr_h2osat",
	2976	& "saturation volume mixing ratio","mol/mol",
	2977	& 3,vmr)
[226]	2978	call wstats(ngrid,"h2o_ice_s",
	2979	& "surface h2o_ice","kg/m2",
[2900]	2980	& 2,qsurf(1,igcm_h2o_ice,iflat))
[758]	2981	call wstats(ngrid,'albedo',
	2982	& 'albedo',
[2900]	2983	& '',2,albedo(1,1,iflat))
[226]	2984	call wstats(ngrid,"mtot",
	2985	& "total mass of water vapor","kg/m2",
	2986	& 2,mtot)
	2987	call wstats(ngrid,"icetot",
	2988	& "total mass of water ice","kg/m2",
	2989	& 2,icetot)
[520]	2990	call wstats(ngrid,"reffice",
	2991	& "Mean reff","m",
	2992	& 2,rave)
[669]	2993	call wstats(ngrid,"Nccntot",
[520]	2994	& "condensation nuclei","Nbr/m2",
[669]	2995	& 2,Nccntot)
	2996	call wstats(ngrid,"Mccntot",
	2997	& "condensation nuclei mass","kg/m2",
	2998	& 2,Mccntot)
[520]	2999	call wstats(ngrid,"rice",
	3000	& "Ice particle size","m",
	3001	& 3,rice)
[226]	3002	if (.not.activice) then
	3003	call wstats(ngrid,"tauTESap",
	3004	& "tau abs 825 cm-1","",
	3005	& 2,tauTES)
[520]	3006	else
[1047]	3007	call wstats(ngrid,'tauTES',
[520]	3008	& 'tau abs 825 cm-1',
	3009	& '',2,taucloudtes)
[226]	3010	endif
	3011
	3012	endif ! of if (water)
[2551]	3013
	3014	if (co2clouds) then
	3015	call wstats(ngrid,"mtotco2",
	3016	& "total mass atm of co2","kg/m2",
	3017	& 2,mtotco2)
	3018	call wstats(ngrid,"icetotco2",
	3019	& "total mass atm of co2 ice","kg/m2",
	3020	& 2,icetotco2)
	3021	call wstats(ngrid,"vaptotco2",
	3022	& "total mass atm of co2 vapor","kg/m2",
[3304]	3023	& 2,vaptotco2)
[2551]	3024	end if
[3203]	3025
	3026
[719]	3027	if (dustbin.ne.0) then
[3203]	3028
[1047]	3029	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
[3203]	3030
[719]	3031	if (doubleq) then
[1047]	3032	call wstats(ngrid,'dqsdust',
[719]	3033	& 'deposited surface dust mass',
	3034	& 'kg.m-2.s-1',2,dqdustsurf)
[1047]	3035	call wstats(ngrid,'dqndust',
[719]	3036	& 'deposited surface dust number',
	3037	& 'number.m-2.s-1',2,dndustsurf)
[1047]	3038	call wstats(ngrid,'reffdust','reffdust',
[719]	3039	& 'm',3,rdust*ref_r0)
[1047]	3040	call wstats(ngrid,'dustq','Dust mass mr',
[756]	3041	& 'kg/kg',3,qdust)
[1047]	3042	call wstats(ngrid,'dustN','Dust number',
[756]	3043	& 'part/kg',3,ndust)
[1974]	3044	if (rdstorm) then
	3045	call wstats(ngrid,'reffstormdust','reffdust',
	3046	& 'm',3,rstormdust*ref_r0)
	3047	call wstats(ngrid,'rdsdustq','Dust mass mr',
	3048	& 'kg/kg',3,rdsqdust)
	3049	call wstats(ngrid,'rdsdustN','Dust number',
	3050	& 'part/kg',3,rdsndust)
	3051	end if
[719]	3052	else
	3053	do iq=1,dustbin
	3054	write(str2(1:2),'(i2.2)') iq
[1047]	3055	call wstats(ngrid,'q'//str2,'mix. ratio',
[719]	3056	& 'kg/kg',3,zq(1,1,iq))
[1047]	3057	call wstats(ngrid,'qsurf'//str2,'qsurf',
[2900]	3058	& 'kg.m-2',2,qsurf(1,iq,iflat))
[719]	3059	end do
	3060	endif ! (doubleq)
[226]	3061
[719]	3062	if (scavenging) then
[1047]	3063	call wstats(ngrid,'ccnq','CCN mass mr',
[756]	3064	& 'kg/kg',3,qccn)
[1047]	3065	call wstats(ngrid,'ccnN','CCN number',
[756]	3066	& 'part/kg',3,nccn)
[719]	3067	endif ! (scavenging)
[3203]	3068
[1464]	3069	endif ! (dustbin.ne.0)
[719]	3070
[2158]	3071	if (photochem) then
[1464]	3072	do iq=1,nq
	3073	if (noms(iq) .ne. "dust_mass" .and.
	3074	$ noms(iq) .ne. "dust_number" .and.
	3075	$ noms(iq) .ne. "ccn_mass" .and.
[1617]	3076	$ noms(iq) .ne. "ccn_number" .and.
	3077	$ noms(iq) .ne. "ccnco2_mass" .and.
[2825]	3078	$ noms(iq) .ne. "ccnco2_number" .and.
	3079	$ noms(iq) .ne. "stormdust_mass" .and.
	3080	$ noms(iq) .ne. "stormdust_number" .and.
	3081	$ noms(iq) .ne. "topdust_mass" .and.
	3082	$ noms(iq) .ne. "topdust_number") then
[1464]	3083	! volume mixing ratio
	3084
[705]	3085	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
[1464]	3086	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
	3087
	3088	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	3089	$ "Volume mixing ratio","mol/mol",3,vmr)
	3090	if ((noms(iq).eq."o")
	3091	$ .or. (noms(iq).eq."co2")
	3092	$ .or. (noms(iq).eq."o3")
	3093	$ .or. (noms(iq).eq."ar")
	3094	$ .or. (noms(iq).eq."o2")
	3095	$ .or. (noms(iq).eq."h2o_vap") ) then
[2932]	3096	call write_output("vmr_"//trim(noms(iq)),
	3097	$ "Volume mixing ratio","mol/mol",vmr(:,:))
[1464]	3098	end if
	3099
	3100	! number density (molecule.cm-3)
	3101
[705]	3102	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	3103	& rho(1:ngrid,1:nlayer)n_avog/
	3104	& (1000*mmol(iq))
[1464]	3105
[2481]	3106	call wstats(ngrid,"num_"//trim(noms(iq)),
[2158]	3107	$ "Number density","cm-3",3,rhopart)
[2932]	3108	call write_output("num_"//trim(noms(iq)),
	3109	$ "Number density","cm-3",rhopart(:,:))
[226]	3110
[1464]	3111	! vertical column (molecule.cm-2)
	3112
	3113	do ig = 1,ngrid
[3203]	3114	colden(ig,iq) = 0.
[1464]	3115	end do
[3203]	3116	do l=1,nlayer
	3117	do ig=1,ngrid
[1464]	3118	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
[3203]	3119	$ *(zplev(ig,l)-zplev(ig,l+1))
	3120	$ 6.022e22/(mmol(iq)g)
	3121	end do
	3122	end do
[1464]	3123
[3203]	3124	call wstats(ngrid,"c_"//trim(noms(iq)),
	3125	$ "column","mol cm-2",2,colden(1,iq))
[2932]	3126	call write_output("c_"//trim(noms(iq)),
	3127	$ "column","mol cm-2",colden(:,iq))
[1464]	3128
	3129	! global mass (g)
[3203]	3130
[1464]	3131	call planetwide_sumval(colden(:,iq)/6.022e23
[1541]	3132	$ mmol(iq)1.e4*cell_area(:),mass(iq))
[1464]	3133
[2932]	3134	call write_output("mass_"//trim(noms(iq)),
	3135	$ "global mass","g",mass(iq))
[1464]	3136
	3137	end if ! of if (noms(iq) .ne. "dust_mass" ...)
	3138	end do ! of do iq=1,nq
[2158]	3139	end if ! of if (photochem)
[1464]	3140
[2563]	3141	IF(lastcall.and.callstats) THEN
[226]	3142	write (,) "Writing stats..."
	3143	call mkstats(ierr)
	3144	ENDIF
	3145
[42]	3146	c (Store EOF for Mars Climate database software)
	3147	IF (calleofdump) THEN
	3148	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	3149	ENDIF
[1212]	3150	#endif
	3151	!endif of ifndef MESOSCALE
[42]	3152
[1236]	3153	#ifdef MESOSCALE
[3203]	3154
	3155	!! see comm_wrf.
[1236]	3156	!! not needed when an array is already in a shared module.
	3157	!! --> example : hfmax_th, zmax_th
[234]	3158
[1590]	3159	CALL allocate_comm_wrf(ngrid,nlayer)
	3160
[1236]	3161	!state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s"
	3162	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
	3163	!state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s"
	3164	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
	3165	!state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2"
[2685]	3166	comm_SWDOWNZ(1:ngrid) = fluxsurf_dn_sw_tot(1:ngrid)
[1236]	3167	!state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" ""
[2415]	3168	comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid)
[1236]	3169	!state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m"
	3170	comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
	3171	!state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2"
[308]	3172	IF (igcm_dust_mass .ne. 0) THEN
[1236]	3173	comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
	3174	ELSE
	3175	comm_QSURFDUST(1:ngrid) = 0.
[308]	3176	ENDIF
[1236]	3177	!state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic"
	3178	comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
	3179	!state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2"
	3180	comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
	3181	!state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
	3182	IF (igcm_h2o_ice .ne. 0) THEN
	3183	comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
	3184	. * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	3185	. * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
	3186	ELSE
	3187	comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
[81]	3188	ENDIF
[1236]	3189	!state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" ""
[1292]	3190	if (activice) then
	3191	comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
	3192	else
	3193	comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
	3194	endif
[1236]	3195	!state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m"
	3196	comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
[3203]	3197
[1242]	3198	!! calculate sensible heat flux in W/m2 for outputs
	3199	!! -- the one computed in vdifc is not the real one
	3200	!! -- vdifc must have been called
	3201	if (.not.callrichsl) then
	3202	sensibFlux(1:ngrid) = zflubid(1:ngrid)
	3203	. - capcal(1:ngrid)*zdtsdif(1:ngrid)
	3204	else
[3203]	3205	sensibFlux(1:ngrid) =
[1242]	3206	& (pplay(1:ngrid,1)/(rpt(1:ngrid,1)))cpp
	3207	& sqrt(pu(1:ngrid,1)pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
	3208	& +(log(1.+0.7wstar(1:ngrid) + 2.3wstar(1:ngrid)2))2)
	3209	& zcdh(1:ngrid)(tsurf(1:ngrid)-zh(1:ngrid,1))
	3210	endif
[3203]	3211
[226]	3212	#else
[528]	3213	#ifndef MESOINI
[42]	3214
[226]	3215	c ==========================================================
	3216	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	3217	c any variable for diagnostic (output with period
	3218	c "ecritphy", set in "run.def")
	3219	c ==========================================================
	3220	c WRITEDIAGFI can ALSO be called from any other subroutines
	3221	c for any variables !!
[2934]	3222	call write_output("emis","Surface emissivity","",
[2900]	3223	& emis(:,iflat))
[2907]	3224	do islope=1,nslope
	3225	write(str2(1:2),'(i2.2)') islope
[2932]	3226	call write_output("emis_slope"//str2,
[2934]	3227	& "Surface emissivity","",emis(:,islope))
[2907]	3228	ENDDO
[2934]	3229	call write_output("zzlay","Midlayer altitude",
	3230	& "m",zzlay(:,:))
[2932]	3231	call write_output("zzlev","Interlayer altitude",
	3232	& "m",zzlev(:,1:nlayer))
	3233	call write_output("pphi","Geopotential","m2s-2",
	3234	& pphi(:,:))
	3235	call write_output("phisfi","Surface geopotential",
	3236	& "m2s-2",phisfi(:))
[3078]	3237	if (grid_type == regular_lonlat) then
	3238	call write_output("area","Mesh area","m2",
	3239	& cell_area_for_lonlat_outputs)
	3240	else ! unstructured grid (e.g. dynamico)
	3241	call write_output("area","Mesh area","m2",cell_area)
	3242	endif
[2932]	3243	call write_output("tsurf","Surface temperature","K",
[2900]	3244	& tsurf(:,iflat))
[2907]	3245	do islope=1,nslope
	3246	write(str2(1:2),'(i2.2)') islope
[2932]	3247	call write_output("tsurf_slope"//str2,
	3248	& "Surface temperature","K",
[2907]	3249	& tsurf(:,islope))
	3250	ENDDO
[2932]	3251	call write_output("ps","surface pressure","Pa",ps(:))
	3252	call write_output("co2ice","co2 ice thickness"
	3253	& ,"kg.m-2",qsurf(:,igcm_co2,iflat))
[2907]	3254	do islope=1,nslope
	3255	write(str2(1:2),'(i2.2)') islope
[2932]	3256	call write_output("co2ice_slope"//str2,"co2 ice thickness"
	3257	& ,"kg.m-2",qsurf(:,igcm_co2,islope))
[2907]	3258	ENDDO
[2934]	3259	call write_output("watercap","Perennial water ice thickness"
[2932]	3260	& ,"kg.m-2",watercap(:,iflat))
[2907]	3261	do islope=1,nslope
	3262	write(str2(1:2),'(i2.2)') islope
[2932]	3263	call write_output("watercap_slope"//str2,
[2934]	3264	& "Perennial water ice thickness"
[2932]	3265	& ,"kg.m-2",watercap(:,islope))
[2907]	3266	ENDDO
[3130]	3267	call write_output("perennial_co2ice",
	3268	& "Perennial co2 ice thickness","kg.m-2",
	3269	& perennial_co2ice(:,iflat))
	3270	do islope=1,nslope
	3271	write(str2(1:2),'(i2.2)') islope
	3272	call write_output("perennial_co2ice_slope"//str2,
[3152]	3273	& "Perennial co2 ice thickness"
[3130]	3274	& ,"kg.m-2",perennial_co2ice(:,islope))
	3275	ENDDO
[2932]	3276	call write_output("temp_layer1","temperature in layer 1",
	3277	& "K",zt(:,1))
	3278	call write_output("temp7","temperature in layer 7",
	3279	& "K",zt(:,7))
	3280	call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2",
[2900]	3281	& fluxsurf_lw(:,iflat))
[2907]	3282	do islope=1,nslope
	3283	write(str2(1:2),'(i2.2)') islope
[2932]	3284	call write_output("fluxsurf_lw_slope"//str2,
	3285	& "fluxsurf_lw","W.m-2",
[2907]	3286	& fluxsurf_lw(:,islope))
	3287	ENDDO
[2932]	3288	call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw",
	3289	& "W.m-2",fluxsurf_dn_sw_tot(:,iflat))
[2907]	3290	do islope=1,nslope
	3291	write(str2(1:2),'(i2.2)') islope
[2932]	3292	call write_output("fluxsurf_dn_sw_slope"//str2,
[2907]	3293	& "fluxsurf_dn_sw",
[2932]	3294	& "W.m-2",fluxsurf_dn_sw_tot(:,islope))
[2907]	3295	ENDDO
[3037]	3296	call write_output("fluxtop_dn_sw","fluxtop_dn_sw",
	3297	& "W.m-2",fluxtop_dn_sw(:,1) + fluxtop_dn_sw(:,2))
[2932]	3298	call write_output("fluxtop_lw","fluxtop_lw","W.m-2",
	3299	& fluxtop_lw(:))
	3300	call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2",
	3301	& fluxtop_up_sw_tot(:))
	3302	call write_output("temp","temperature","K",zt(:,:))
	3303	call write_output("Sols","Time","sols",zday)
	3304	call write_output("Ls","Solar longitude","deg",
	3305	& zls*180./pi)
	3306	call write_output("u","Zonal wind","m.s-1",zu(:,:))
	3307	call write_output("v","Meridional wind","m.s-1",zv(:,:))
	3308	call write_output("w","Vertical wind","m.s-1",pw(:,:))
	3309	call write_output("rho","density","kg.m-3",rho(:,:))
	3310	call write_output("pressure","Pressure","Pa",zplay(:,:))
[2934]	3311	call write_output("zplev","Interlayer pressure","Pa",
	3312	& zplev(:,1:nlayer))
[2932]	3313	call write_output('sw_htrt','sw heat. rate',
	3314	& 'K/s',zdtsw(:,:))
	3315	call write_output('lw_htrt','lw heat. rate',
	3316	& 'K/s',zdtlw(:,:))
	3317	call write_output("local_time","Local time",
	3318	& 'sol',local_time(:))
[3147]	3319	if (water) then
[520]	3320	if (.not.activice) then
[2932]	3321	CALL write_output('tauTESap',
[520]	3322	& 'tau abs 825 cm-1',
[2932]	3323	& '',tauTES(:))
[520]	3324	else
[2932]	3325	CALL write_output('tauTES',
[520]	3326	& 'tau abs 825 cm-1',
[2932]	3327	& '',taucloudtes(:))
[520]	3328	endif
[3147]	3329	endif ! of if (water)
[528]	3330	#else
[2362]	3331	!!! this is to ensure correct initialisation of mesoscale model
[2932]	3332	call write_output("tsurf","Surface temperature","K",
[2900]	3333	& tsurf(:,iflat))
[2932]	3334	call write_output("ps","surface pressure","Pa",ps(:))
	3335	call write_output("co2ice","co2 ice thickness","kg.m-2",
[2900]	3336	& qsurf(:,igcm_co2,iflat))
[2932]	3337	call write_output("temp","temperature","K",zt(:,:))
	3338	call write_output("u","Zonal wind","m.s-1",zu(:,:))
	3339	call write_output("v","Meridional wind","m.s-1",zv(:,:))
[2934]	3340	call write_output("emis","Surface emissivity","",
[2932]	3341	& emis(:,iflat))
	3342	call write_output("tsoil","Soil temperature",
	3343	& "K",tsoil(:,:,iflat))
	3344	call write_output("inertiedat","Soil inertia",
	3345	& "K",inertiedat(:,:))
[299]	3346	#endif
[42]	3347
[226]	3348	c ----------------------------------------------------------
	3349	c Outputs of the CO2 cycle
	3350	c ----------------------------------------------------------
[2823]	3351	if (igcm_co2.ne.0) then
[2932]	3352	call write_output("co2","co2 mass mixing ratio",
	3353	& "kg.kg-1",zq(:,:,igcm_co2))
[2362]	3354
	3355	if (co2clouds) then
[2932]	3356	call write_output('ccnqco2','CCNco2 mmr',
	3357	& 'kg.kg-1',zq(:,:,igcm_ccnco2_mass))
[1660]	3358
[2932]	3359	call write_output('ccnNco2','CCNco2 number',
	3360	& 'part.kg-1',zq(:,:,igcm_ccnco2_number))
[2362]	3361
[2932]	3362	call write_output('co2_ice','co2_ice mmr in atm',
	3363	& 'kg.kg-1',zq(:,:,igcm_co2_ice))
[2362]	3364
[2932]	3365	call write_output("mtotco2","total mass atm of co2",
	3366	& "kg.m-2",mtotco2(:))
	3367	call write_output("icetotco2","total mass atm of co2 ice",
	3368	& "kg.m-2", icetotco2(:))
	3369	call write_output("vaptotco2","total mass atm of co2 "//
	3370	& "vapor","kg.m-2", vaptotco2(:))
[2589]	3371	if (co2useh2o) then
[2932]	3372	call write_output('ccnqco2_h2o_m_ice',
[2589]	3373	& 'CCNco2_h2o_mass_ice mmr',
[2932]	3374	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice))
[2589]	3375
[2932]	3376	call write_output('ccnqco2_h2o_m_ccn',
[2589]	3377	& 'CCNco2_h2o_mass_ccn mmr',
[2932]	3378	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn))
[2589]	3379
[2932]	3380	call write_output('ccnNco2_h2o','CCNco2_h2o number',
	3381	& 'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number))
[2589]	3382	end if
	3383
	3384	if (meteo_flux) then
[2932]	3385	call write_output('ccnqco2_meteor','CCNco2_meteor mmr',
	3386	& 'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass))
[2589]	3387
[2932]	3388	call write_output('ccnNco2_meteor','CCNco2_meteor number',
	3389	& 'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number))
[2589]	3390	end if
	3391
[2362]	3392	end if ! of if (co2clouds)
[2823]	3393	end if ! of if (igcm_co2.ne.0)
[2362]	3394
	3395	! Output He tracer, if there is one
[2823]	3396	if (igcm_he.ne.0) then
[2932]	3397	call write_output("he","helium mass mixing ratio",
	3398	& "kg/kg",zq(:,:,igcm_he))
[2362]	3399	vmr = zq(1:ngrid,1:nlayer,igcm_he)
	3400	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
[2932]	3401	call write_output('vmr_he','helium vol. mixing ratio',
	3402	& 'mol/mol',vmr(:,:))
[2362]	3403	end if
	3404
[226]	3405	c ----------------------------------------------------------
	3406	c Outputs of the water cycle
	3407	c ----------------------------------------------------------
[2362]	3408	if (water) then
[299]	3409	#ifdef MESOINI
	3410	!!!! waterice = q01, voir readmeteo.F90
[2932]	3411	call write_output('q01',noms(igcm_h2o_ice),
	3412	& 'kg/kg',
	3413	& zq(:,:,igcm_h2o_ice))
[299]	3414	!!!! watervapor = q02, voir readmeteo.F90
[2932]	3415	call write_output('q02',noms(igcm_h2o_vap),
	3416	& 'kg/kg',
	3417	& zq(:,:,igcm_h2o_vap))
[299]	3418	!!!! surface waterice qsurf02 (voir readmeteo)
[2932]	3419	call write_output('qsurf02','surface tracer',
	3420	& 'kg.m-2',
	3421	& qsurf(:,igcm_h2o_ice,iflat))
[2907]	3422	do islope=1,nslope
	3423	write(str2(1:2),'(i2.2)') islope
[2932]	3424	call write_output('qsurf02_slope'//str2,
	3425	& 'surface tracer','kg.m-2',
	3426	& qsurf(:,igcm_h2o_ice,islope))
[2907]	3427	ENDDO
[299]	3428	#endif
[2932]	3429	call write_output('mtot',
[2362]	3430	& 'total mass of water vapor',
[2932]	3431	& 'kg/m2',mtot(:))
	3432	call write_output('icetot',
[2362]	3433	& 'total mass of water ice',
[2932]	3434	& 'kg/m2',icetot(:))
[2362]	3435	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	3436	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
[2932]	3437	call write_output('vmr_h2oice','h2o ice vmr',
	3438	& 'mol/mol',vmr(:,:))
[2362]	3439	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
	3440	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
[2932]	3441	call write_output('vmr_h2ovap','h2o vap vmr',
	3442	& 'mol/mol',vmr(:,:))
	3443	call write_output('reffice',
[2362]	3444	& 'Mean reff',
[2932]	3445	& 'm',rave(:))
	3446	call write_output('h2o_ice','h2o_ice','kg/kg',
	3447	& zq(:,:,igcm_h2o_ice))
	3448	call write_output('h2o_vap','h2o_vap','kg/kg',
	3449	& zq(:,:,igcm_h2o_vap))
[1922]	3450
[2312]	3451	if (hdo) then
	3452	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
	3453	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
[2932]	3454	CALL write_output('vmr_hdoice','hdo ice vmr',
	3455	& 'mol/mol',vmr(:,:))
[2312]	3456	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
	3457	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
[2932]	3458	CALL write_output('vmr_hdovap','hdo vap vmr',
	3459	& 'mol/mol',vmr(:,:))
	3460	call write_output('hdo_ice','hdo_ice','kg/kg',
	3461	& zq(:,:,igcm_hdo_ice))
	3462	call write_output('hdo_vap','hdo_vap','kg/kg',
	3463	& zq(:,:,igcm_hdo_vap))
[1922]	3464
[2932]	3465	CALL write_output('mtotD',
[2312]	3466	& 'total mass of HDO vapor',
[2932]	3467	& 'kg/m2',mtotD(:))
	3468	CALL write_output('icetotD',
[2312]	3469	& 'total mass of HDO ice',
[2932]	3470	& 'kg/m2',icetotD(:))
[2312]	3471
	3472	C Calculation of the D/H ratio
	3473	do l=1,nlayer
	3474	do ig=1,ngrid
[2362]	3475	if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
[2312]	3476	DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
	3477	& zq(ig,l,igcm_h2o_vap) )1./(2.155.76e-6)
	3478	else
	3479	DoH_vap(ig,l) = 0.
	3480	endif
	3481	enddo
	3482	enddo
	3483
	3484	do l=1,nlayer
	3485	do ig=1,ngrid
[2324]	3486	if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then
[2312]	3487	DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
	3488	& zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
	3489	else
	3490	DoH_ice(ig,l) = 0.
	3491	endif
	3492	enddo
	3493	enddo
	3494
[2932]	3495	CALL write_output('DoH_vap',
[2312]	3496	& 'D/H ratio in vapor',
[2932]	3497	& ' ',DoH_vap(:,:))
	3498	CALL write_output('DoH_ice',
[2312]	3499	& 'D/H ratio in ice',
[2932]	3500	& '',DoH_ice(:,:))
[2312]	3501
	3502	endif !hdo
	3503
[1711]	3504	!A. Pottier
[2934]	3505	! CALL write_output('rmoym',
	3506	! & 'alternative reffice',
	3507	! & 'm',rave2(:))
	3508	call write_output('h2o_saturation',
	3509	& 'h2o vap saturation ratio','',satu(:,:))
[1130]	3510	if (scavenging) then
[2932]	3511	CALL write_output("Nccntot",
[669]	3512	& "condensation nuclei","Nbr/m2",
[2932]	3513	& Nccntot(:))
	3514	CALL write_output("Mccntot",
[719]	3515	& "mass condensation nuclei","kg/m2",
[2932]	3516	& Mccntot(:))
[1130]	3517	endif
[2934]	3518	call write_output('rice','Water ice particle size',
[2932]	3519	& 'm',rice(:,:))
	3520	call write_output('h2o_ice_s',
[226]	3521	& 'surface h2o_ice',
[2932]	3522	& 'kg.m-2',qsurf(:,igcm_h2o_ice,iflat))
[2907]	3523	do islope=1,nslope
	3524	write(str2(1:2),'(i2.2)') islope
[2932]	3525	call write_output('h2o_ice_s_slope'//str2,
[2907]	3526	& 'surface h2o_ice',
[2932]	3527	& 'kg.m-2',qsurf(:,igcm_h2o_ice,islope))
[2907]	3528	ENDDO
[2312]	3529	if (hdo) then
[2932]	3530	call write_output('hdo_ice_s',
[2312]	3531	& 'surface hdo_ice',
[2932]	3532	& 'kg.m-2',qsurf(:,igcm_hdo_ice,iflat))
[2907]	3533	do islope=1,nslope
	3534	write(str2(1:2),'(i2.2)') islope
[2932]	3535	call write_output('hdo_ice_s_slope'//str2,
[2907]	3536	& 'surface hdo_ice',
[2932]	3537	& 'kg.m-2',qsurf(:,igcm_hdo_ice,islope))
[2907]	3538	ENDDO
[2312]	3539
	3540	do ig=1,ngrid
[2900]	3541	if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
	3542	DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
	3543	& qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
[2312]	3544	else
	3545	DoH_surf(ig) = 0.
	3546	endif
	3547	enddo
	3548
[2932]	3549	call write_output('DoH_surf',
[2312]	3550	& 'surface D/H',
[2932]	3551	& '',DoH_surf(:))
[2312]	3552	endif ! hdo
	3553
[2932]	3554	CALL write_output('albedo',
[833]	3555	& 'albedo',
[2932]	3556	& '',albedo(:,1,iflat))
[2907]	3557	do islope=1,nslope
	3558	write(str2(1:2),'(i2.2)') islope
[2932]	3559	CALL write_output('albedo_slope'//str2,
[2907]	3560	& 'albedo',
[2932]	3561	& '',albedo(:,1,islope))
[2907]	3562	ENDDO
[3230]	3563	if (surfaceice_tifeedback.or.poreice_tifeedback) then
[2932]	3564	call write_output("soiltemp",
[833]	3565	& "Soil temperature","K",
[2932]	3566	& tsoil(:,:,iflat))
	3567	call write_output('soilti',
[833]	3568	& 'Soil Thermal Inertia',
[2942]	3569	& 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
[3230]	3570
[2907]	3571	do islope=1,nslope
	3572	write(str2(1:2),'(i2.2)') islope
[2932]	3573	call write_output('soilti_slope'//str2,
[2907]	3574	& 'Soil Thermal Inertia',
[2942]	3575	& 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope))
[2907]	3576	ENDDO
[833]	3577	endif
[1711]	3578	!A. Pottier
[2362]	3579	if (CLFvarying) then !AP14 nebulosity
[2932]	3580	call write_output('totcloudfrac',
[2362]	3581	& 'Total cloud fraction',
[2932]	3582	& ' ',totcloudfrac(:))
[2362]	3583	end if !clf varying
	3584	end if !(water)
[1711]	3585
[226]	3586	c ----------------------------------------------------------
	3587	c Outputs of the dust cycle
	3588	c ----------------------------------------------------------
	3589
[2932]	3590	call write_output('tau_pref_scenario',
[2415]	3591	& 'Prescribed visible dust optical depth at 610Pa',
[2932]	3592	& 'NU',tau_pref_scenario(:))
[226]	3593
[2932]	3594	call write_output('tau_pref_gcm',
[2415]	3595	& 'Visible dust optical depth at 610Pa in the GCM',
[2932]	3596	& 'NU',tau_pref_gcm(:))
[3203]	3597
[2643]	3598	if (reff_driven_IRtoVIS_scenario) then
[2932]	3599	call write_output('IRtoVIScoef',
[2934]	3600	& 'Conversion coeff for dust tau from abs9.3um to ext0.67um',
[2932]	3601	& '/',IRtoVIScoef(:))
[2643]	3602	endif
[2415]	3603
[2823]	3604	if (dustbin.ne.0) then
[1264]	3605
	3606	#ifndef MESOINI
[226]	3607	if (doubleq) then
[2932]	3608	call write_output('dqsdust',
[719]	3609	& 'deposited surface dust mass',
[2932]	3610	& 'kg.m-2.s-1',dqdustsurf(:))
	3611	call write_output('dqndust',
[719]	3612	& 'deposited surface dust number',
[2932]	3613	& 'number.m-2.s-1',dndustsurf(:))
	3614	call write_output('reffdust','reffdust',
	3615	& 'm',rdust(:,:)*ref_r0)
	3616	call write_output('dustq','Dust mass mr',
	3617	& 'kg/kg',qdust(:,:))
	3618	call write_output('dustN','Dust number',
	3619	& 'part/kg',ndust(:,:))
[2362]	3620
[3203]	3621	select case (trim(dustiropacity))
[2247]	3622	case ("tes")
[2934]	3623	call write_output('dsodust_TES',
[2247]	3624	& 'density scaled extinction opacity of std dust at 9.3um(TES)',
[2932]	3625	& 'm2.kg-1',dsodust(:,:))
[2934]	3626	call write_output('dso_TES',
[2247]	3627	& 'density scaled extinction opacity of all dust at 9.3um(TES)',
[2932]	3628	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
[2247]	3629	case ("mcs")
[2932]	3630	call write_output('dsodust',
[2247]	3631	& 'density scaled extinction opacity of std dust at 21.6um(MCS)',
[2932]	3632	& 'm2.kg-1',dsodust(:,:))
	3633	call write_output('dso',
[2247]	3634	& 'density scaled extinction opacity of all dust at 21.6um(MCS)',
[2932]	3635	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
[2247]	3636	end select
[2246]	3637	else ! (doubleq=.false.)
[226]	3638	do iq=1,dustbin
	3639	write(str2(1:2),'(i2.2)') iq
[2932]	3640	call write_output('q'//str2,'mix. ratio',
	3641	& 'kg/kg',zq(:,:,iq))
	3642	call write_output('qsurf'//str2,'qsurf',
	3643	& 'kg.m-2',qsurf(:,iq,iflat))
[2907]	3644	do islope=1,nslope
	3645	write(str2(1:2),'(i2.2)') islope
[2932]	3646	call write_output('qsurf_slope'//str2,'qsurf',
	3647	& 'kg.m-2',qsurf(:,iq,islope))
[2907]	3648	ENDDO
[226]	3649	end do
	3650	endif ! (doubleq)
[358]	3651
[1974]	3652	if (rdstorm) then ! writediagfi tendencies stormdust tracers
[2932]	3653	call write_output('reffstormdust','reffstormdust',
	3654	& 'm',rstormdust(:,:)*ref_r0)
	3655	call write_output('mstormdtot',
[1974]	3656	& 'total mass of stormdust only',
[2932]	3657	& 'kg.m-2',mstormdtot(:))
	3658	call write_output('mdusttot',
[1974]	3659	& 'total mass of dust only',
[2932]	3660	& 'kg.m-2',mdusttot(:))
	3661	call write_output('rdsdqsdust',
[1974]	3662	& 'deposited surface stormdust mass',
[2932]	3663	& 'kg.m-2.s-1',rdsdqdustsurf(:))
	3664	call write_output('rdsdustq','storm Dust mass mr',
	3665	& 'kg/kg',rdsqdust(:,:))
	3666	call write_output('rdsdustqmodel','storm Dust massmr',
	3667	& 'kg/kg',pq(:,:,igcm_stormdust_mass))
	3668	call write_output('rdsdustN','storm Dust number',
	3669	& 'part/kg',rdsndust(:,:))
	3670	call write_output("stormfract",
[1974]	3671	& "fraction of the mesh, with stormdust","none",
[2932]	3672	& totstormfract(:))
[2934]	3673	! call write_output('qsurf',
	3674	! & 'stormdust injection',
	3675	! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat))
	3676	! do islope=1,nslope
	3677	! write(str2(1:2),'(i2.2)') islope
	3678	! call write_output('qsurf_slope'//str2,
	3679	! & 'stormdust injection',
	3680	! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,islope))
	3681	! ENDDO
	3682	! call write_output('pdqsurf',
	3683	! & 'tendancy stormdust mass at surface',
	3684	! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat))
	3685	! do islope=1,nslope
	3686	! write(str2(1:2),'(i2.2)') islope
	3687	! call write_output('pdqsurf_slope'//str2,
	3688	! & 'tendancy stormdust mass at surface',
	3689	! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope))
	3690	! ENDDO
	3691	call write_output('wspeed_stormdust',
	3692	& 'vertical velocity of stormdust',
	3693	& 'm/s',wspeed(:,:))
	3694	call write_output('zdqsed_dust_mass'
	3695	& ,'sedimentation tendency of background dust mmr'
	3696	& ,'kg/kg.s-1',
[1974]	3697	& zdqsed(:,:,igcm_dust_mass))
[2934]	3698	call write_output('zdqssed_dust_mass'
	3699	& ,'sedimentation tendency of background dust on surface'
	3700	& ,'kg.m-2.s-1',
[1974]	3701	& zdqssed(:,igcm_dust_mass))
[2934]	3702	call write_output('zdqsed_stormdust_mass'
	3703	& ,'sedimentation tendency of stormdust dust mmr'
	3704	& ,'kg/kg.s-1',
[1974]	3705	& zdqsed(:,:,igcm_stormdust_mass))
[2934]	3706	call write_output('zdqsed_dust_number'
	3707	& ,'sedimentation tendency of background dust number'
	3708	& ,'nbr/kg.s-1',
[2931]	3709	& zdqsed(:,:,igcm_dust_number))
[2932]	3710	call write_output('rdust','rdust',
	3711	& 'm',rdust(:,:))
	3712	call write_output('rstormdust','rstormdust',
	3713	& 'm',rstormdust(:,:))
[2362]	3714
[3203]	3715	select case (trim(dustiropacity))
[2247]	3716	case ("tes")
[2934]	3717	call write_output('dsords_TES',
[2247]	3718	& 'density scaled extinction opacity of stormdust at 9.3um(TES)',
[2932]	3719	& 'm2.kg-1',dsords(:,:))
[2247]	3720	case ("mcs")
[2932]	3721	call write_output('dsords',
[2247]	3722	& 'density scaled extinction opacity of stormdust at 21.6um(MCS)',
[2932]	3723	& 'm2.kg-1',dsords(:,:))
[2247]	3724	end select
[1974]	3725	endif ! (rdstorm)
	3726
[2628]	3727	if (topflows) then
[2934]	3728	call write_output('refftopdust',
	3729	& 'Topdust dust effective radius',
	3730	& 'm',rtopdust(:,:)*ref_r0)
[2932]	3731	call write_output('topdustq','top Dust mass mr',
	3732	& 'kg/kg',pq(:,:,igcm_topdust_mass))
	3733	call write_output('topdustN','top Dust number',
	3734	& 'part/kg',pq(:,:,igcm_topdust_number))
[3203]	3735	select case (trim(dustiropacity))
[2247]	3736	case ("tes")
[2934]	3737	call write_output('dsotop_TES',
[2247]	3738	& 'density scaled extinction opacity of topdust at 9.3um(TES)',
[2932]	3739	& 'm2.kg-1',dsotop(:,:))
[2247]	3740	case ("mcs")
[2932]	3741	call write_output('dsotop',
[2247]	3742	& 'density scaled extinction opacity of topdust at 21.6um(MCS)',
[2932]	3743	& 'm2.kg-1',dsotop(:,:))
[2247]	3744	end select
[2628]	3745	endif ! (topflows)
[2362]	3746
[2417]	3747	if (dustscaling_mode==2) then
[2932]	3748	call write_output("dust_rad_adjust",
[2417]	3749	& "radiative adjustment coefficient for dust",
[2932]	3750	& "",dust_rad_adjust(:))
[2417]	3751	endif
	3752
[2934]	3753	! if (scavenging) then ! these outputs should be in the scavenging routine
	3754	! call write_output('ccnq','CCN mass mr',
	3755	! & 'kg/kg',qccn(:,:))
	3756	! call write_output('ccnN','CCN number',
	3757	! & 'part/kg',nccn(:,:))
	3758	! call write_output('surfccnq','Surf nuclei mass mr',
	3759	! & 'kg.m-2',qsurf(:,igcm_ccn_mass,iflat))
	3760	! do islope=1,nslope
	3761	! write(str2(1:2),'(i2.2)') islope
	3762	! call write_output('surfccnq_slope'//str2,
	3763	! & 'Surf nuclei mass mr',
	3764	! & 'kg.m-2',qsurf(:,igcm_ccn_mass,islope))
	3765	! ENDDO
	3766	! call write_output('surfccnN','Surf nuclei number',
	3767	! & 'kg.m-2',qsurf(:,igcm_ccn_number,iflat))
	3768	! do islope=1,nslope
	3769	! write(str2(1:2),'(i2.2)') islope
	3770	! call write_output('surfccnN_slope'//str2,
	3771	! & 'Surf nuclei number',
	3772	! & 'kg.m-2',qsurf(:,igcm_ccn_number,islope))
	3773	! ENDDO
	3774	! endif ! (scavenging)
[358]	3775
[1264]	3776	#else
	3777	! !!! to initialize mesoscale we need scaled variables
	3778	! !!! because this must correspond to starting point for tracers
[2932]	3779	! call write_output('dustq','Dust mass mr',
[1264]	3780	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
[2932]	3781	! call write_output('dustN','Dust number',
[1264]	3782	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
[2932]	3783	! call write_output('ccn','Nuclei mass mr',
[1264]	3784	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
[2932]	3785	! call write_output('ccnN','Nuclei number',
[1264]	3786	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
	3787	if (freedust) then
[2932]	3788	call write_output('dustq','Dust mass mr',
	3789	& 'kg/kg',qdust)
	3790	call write_output('dustN','Dust number',
	3791	& 'part/kg',ndust)
	3792	call write_output('ccn','CCN mass mr',
	3793	& 'kg/kg',qccn)
	3794	call write_output('ccnN','CCN number',
	3795	& 'part/kg',nccn)
[1264]	3796	else
[2932]	3797	call write_output('dustq','Dust mass mr',
	3798	& 'kg/kg',pq(:,:,igcm_dust_mass))
	3799	call write_output('dustN','Dust number',
	3800	& 'part/kg',pq(:,:,igcm_dust_number))
	3801	call write_output('ccn','Nuclei mass mr',
	3802	& 'kg/kg',pq(:,:,igcm_ccn_mass))
	3803	call write_output('ccnN','Nuclei number',
	3804	& 'part/kg',pq(:,:,igcm_ccn_number))
[1264]	3805	endif
	3806	#endif
	3807
[2823]	3808	end if ! (dustbin.ne.0)
[226]	3809
[2149]	3810	c ----------------------------------------------------------
[705]	3811	c Thermospheric outputs
[267]	3812	c ----------------------------------------------------------
[2362]	3813	if(callthermos) then
[705]	3814
[2932]	3815	call write_output("quv","UV heating","K/s",
	3816	$ zdteuv(:,:))
	3817	call write_output("cond","Thermal conduction","K/s",
	3818	$ zdtconduc(:,:))
[705]	3819
[2467]	3820	!H, H2 and D escape fluxes
	3821
[2932]	3822	call write_output("PhiH","H escape flux","s-1",
	3823	$ PhiEscH)
	3824	call write_output("PhiH2","H2 escape flux","s-1",
	3825	$ PhiEscH2)
	3826	call write_output("PhiD","D escape flux","s-1",
	3827	$ PhiEscD)
[2467]	3828
[705]	3829	endif !(callthermos)
	3830
[2932]	3831	call write_output("q15um","15 um cooling","K/s",
	3832	$ zdtnlte(:,:))
	3833	call write_output("qnir","NIR heating","K/s",
	3834	$ zdtnirco2(:,:))
[2362]	3835
[705]	3836	c ----------------------------------------------------------
	3837	c ----------------------------------------------------------
[267]	3838	c PBL OUTPUS
	3839	c ----------------------------------------------------------
	3840	c ----------------------------------------------------------
	3841
	3842	c ----------------------------------------------------------
[226]	3843	c Outputs of thermals
	3844	c ----------------------------------------------------------
[2362]	3845	if (calltherm) then
[2932]	3846	call write_output('lmax_th',
[2934]	3847	& 'index of vertical extension of thermals',
	3848	& 'grid level',lmax_th_out(:))
[2932]	3849	call write_output('zmax_th',
[2934]	3850	& 'vertical extension of thermals','m',
[2932]	3851	& zmax_th(:))
	3852	call write_output('hfmax_th',
[2934]	3853	& 'maximum heat flux in thermals','K.m/s',
[2932]	3854	& hfmax_th(:))
	3855	call write_output('wstar',
[2934]	3856	& 'maximum thermals vertical velocity','m/s',
[2932]	3857	& wstar(:))
[2362]	3858	end if
[226]	3859
	3860	c ----------------------------------------------------------
[267]	3861	c ----------------------------------------------------------
	3862	c END OF PBL OUTPUS
	3863	c ----------------------------------------------------------
	3864	c ----------------------------------------------------------
	3865
	3866	c ----------------------------------------------------------
[226]	3867	c Output in netcdf file "diagsoil.nc" for subterranean
	3868	c variables (output every "ecritphy", as for writediagfi)
	3869	c ----------------------------------------------------------
	3870	! Write soil temperature
[2932]	3871	call write_output("soiltemp","Soil temperature","K",
	3872	& tsoil(:,:,iflat))
[2907]	3873	do islope=1,nslope
	3874	write(str2(1:2),'(i2.2)') islope
[2932]	3875	call write_output("soiltemp_slope"//str2,
[2907]	3876	& "Soil temperature","K",
[2932]	3877	& tsoil(:,:,islope))
[2907]	3878	ENDDO
[226]	3879
[2931]	3880	!PREVIOUSLY IN 1D ONLY
[2932]	3881	call write_output("dtrad","rad. heat. rate",
	3882	& "K.s-1",dtrad(:,:))
[544]	3883
[1974]	3884	if (rdstorm) then
[2932]	3885	call write_output('aerosol_dust','opacity of env. dust',''
	3886	& ,aerosol(:,:,iaer_dust_doubleq))
	3887	call write_output('aerosol_stormdust',
[1974]	3888	& 'opacity of storm dust',''
[2932]	3889	& ,aerosol(:,:,iaer_stormdust_doubleq))
[2934]	3890	call write_output('dqsdifdustq',
	3891	&'tendency due to vertical diffusion of background dust on surface'
	3892	& ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat))
[2907]	3893	do islope=1,nslope
	3894	write(str2(1:2),'(i2.2)') islope
[2932]	3895	call write_output('dqsdifdustq_slope'//str2,
[2934]	3896	&'tendency due to vertical diffusion of background dust on surface'
	3897	& ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope))
[2907]	3898	ENDDO
[2934]	3899	call write_output('dqsdifrdsq',
	3900	& 'tendency due to vertical diffusion of stormdust on surface',
[2932]	3901	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat))
[2907]	3902	do islope=1,nslope
	3903	write(str2(1:2),'(i2.2)') islope
[2932]	3904	call write_output('dqsdifrdsq_slope'//str2,
[2934]	3905	& 'tendency due to vertical diffusion of stormdust on surface',
[2932]	3906	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope))
[2907]	3907	ENDDO
[2931]	3908	endif !(rdstorm)
[1974]	3909
[2931]	3910	if(water) then
	3911	if (.not.scavenging) then
[2932]	3912	call write_output('zdqcloud_ice','cloud ice',
	3913	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice))
	3914	call write_output('zdqcloud_vap','cloud vap',
	3915	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap))
[3000]	3916	call write_output('zdqcloud','cloud',
[2932]	3917	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)
	3918	& +zdqcloud(:,:,igcm_h2o_vap))
[2312]	3919	IF (hdo) THEN
[2932]	3920	call write_output('zdqcloud_iceD','cloud ice hdo',
	3921	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice))
	3922	call write_output('zdqcloud_vapD','cloud vap hdo',
	3923	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap))
[2312]	3924	ENDIF ! hdo
[2931]	3925	endif !not.scavenging
[2963]	3926
	3927	! Output needed by the PEM
	3928	DO ig = 1,ngrid
	3929	ztmp1 =(1/m_co2 - 1/m_noco2)
	3930	ztmp2=1/m_noco2
	3931	pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2)
	3932	& * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)1.e-3)ps(ig)
	3933
	3934	DO islope = 1,nslope
[3203]	3935	! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
[2963]	3936	rhowater_surf_sat(ig,islope) =
[2965]	3937	& exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o)
[2963]	3938	& / tsurf(ig,islope)
[2965]	3939	& * mmol(igcm_h2o_vap)/(mugaz*r)
[3203]	3940
[2963]	3941	if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then
[3203]	3942	! we consider to be at saturation above 1.e-4 kg.m-2
[2965]	3943	rhowater_surf(ig,islope) = rhowater_surf_sat(ig,islope)
[2963]	3944	else
[3203]	3945	! otherwise, use vapor partial pressure
[2963]	3946	rhowater_surf(ig,islope) = pvap_surf(ig)
	3947	& / tsurf(ig,islope)
[2965]	3948	& * mmol(igcm_h2o_vap)/(mugaz*r)
[2963]	3949	endif
	3950	DO isoil = 1,nsoilmx
	3951	rhowater_soil(ig,isoil,islope) =
[2965]	3952	& exp(beta_clap_h2o/tsoil(ig,isoil,islope)+alpha_clap_h2o)
[2963]	3953	& / tsoil(ig,isoil,islope)
[2965]	3954	& * mmol(igcm_h2o_vap)/(mugaz*r)
[2963]	3955	ENDDO
	3956	ENDDO
	3957	ENDDO
	3958
[3122]	3959	CALL write_output("waterdensity_soil",
[3106]	3960	& "rhowater_soil",'kg.m-3',
	3961	& rhowater_soil(:,:,iflat))
[3122]	3962	CALL write_output("waterdensity_surface",
[3106]	3963	& "rhowater_surface",'kg.m-3',
	3964	& rhowater_surf(:,iflat))
[2963]	3965	DO islope = 1,nslope
	3966	write(str2(1:2),'(i2.2)') islope
[3122]	3967	CALL write_output("waterdensity_soil_slope"//str2,
[2965]	3968	& "rhowater_soil_slope"//str2,'kg.m-3',
[2963]	3969	& rhowater_soil(:,:,islope))
[3612]	3970	CALL write_output("waterdensity_surface_slope"//str2,
[2965]	3971	& "rhowater_surface"//str2,'kg.m-3',
[2963]	3972	& rhowater_surf(:,islope))
	3973	ENDDO
	3974
[2965]	3975	CALL write_output("h2o_layer1","h2o mass mr in the first layer",
	3976	& 'kg/kg',zq(:,1,igcm_h2o_vap))
	3977	CALL write_output("co2_layer1","co2 mass mr in the first layer",
	3978	& 'kg/kg',zq(:,1,igcm_co2))
[3027]	3979	ENDIF ! of IF (water)
[2963]	3980
[2931]	3981	!PREVIOUSLY IN 1D ONLY
[2383]	3982
[2931]	3983	c ==========================================================
	3984	c END OF WRITEDIAGFI
	3985	c ==========================================================
	3986	#endif
	3987	! of ifdef MESOSCALE
[1711]	3988
[2931]	3989	c ELSE ! if(ngrid.eq.1)
[358]	3990
[2931]	3991	c#ifndef MESOSCALE
	3992	c write(*,
	3993	c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
	3994	c & zls*180./pi,odpref,tau_pref_scenario
	3995	c#endif
[226]	3996
[2931]	3997	c END IF ! if(ngrid.ne.1)
[2887]	3998
[2374]	3999	! test for co2 conservation with co2 microphysics
	4000	if (igcm_co2_ice.ne.0) then
	4001	co2totB = 0. ! added by C.M.
	4002	do ig=1,ngrid
	4003	do l=1,nlayer
	4004	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
	4005	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
	4006	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
	4007	enddo
[2900]	4008	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
[2362]	4009	enddo
[2599]	4010	else
	4011	co2totB = 0. ! added by C.M.
	4012	do ig=1,ngrid
	4013	do l=1,nlayer
	4014	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
	4015	& (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
	4016	enddo
[2900]	4017	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
[2599]	4018	enddo
	4019	endif ! of if (igcm_co2_ice.ne.0)
[2660]	4020	co2conservation = (co2totA-co2totB)/co2totA
[2932]	4021	call write_output( 'co2conservation',
[2374]	4022	& 'Total CO2 mass conservation in physic',
[2934]	4023	& 'kg', co2conservation)
[2223]	4024	! XIOS outputs
[3203]	4025	#ifdef CPP_XIOS
[2223]	4026	! Send fields to XIOS: (NB these fields must also be defined as
	4027	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[2545]	4028
	4029	CALL send_xios_field("controle",tab_cntrl_mod,1)
	4030
	4031	CALL send_xios_field("ap",ap,1)
	4032	CALL send_xios_field("bp",bp,1)
	4033	CALL send_xios_field("aps",aps,1)
	4034	CALL send_xios_field("bps",bps,1)
	4035
[2223]	4036	if (lastcall.and.is_omp_master) then
	4037	write(,) "physiq lastcall: call xios_context_finalize"
	4038	call xios_context_finalize
	4039	endif
	4040	#endif
	4041
[2551]	4042	if (check_physics_outputs) then
	4043	! Check the validity of updated fields at the end of the physics step
[2570]	4044	call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq)
[2551]	4045	endif
[2223]	4046
[42]	4047	icount=icount+1
[1549]	4048
	4049	END SUBROUTINE physiq
	4050
	4051	END MODULE physiq_mod

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

Download in other formats:

Original Format