Context Navigation

physiq_mod.F @ 3242

Last change on this file since 3242 was 3230, checked in by llange, 17 months ago

Mars PCM

Move soil_tifeedback into a module waterice_tifeedback_mod.F90
The flag to call tifeedback has been changed from "tifeedback" to surfaceice_tifeedback for clarity
Add the possibility to change the thermal inertia while pore ices is forming in the soil: to do so, use the flag poreice_tifeedback. The computation is done in waterice_tifeedback_mod.F90.

For now, surfaceice_tifeedback and poreice_tifeedback can not be called together, we might think about how to merge later.
LL

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

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

Download in other formats:

Original Format