Context Navigation

physiq_mod.F @ 2917

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

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

Download in other formats:

Original Format