Context Navigation

physiq_mod.F @ 2814

Last change on this file since 2814 was 2808, checked in by emillour, 3 years ago
Mars GCM: Bug fix in jthermcalc_e107 (wrong indexes used in some computations involving NO2 and H2). AM
File size: 163.4 KB

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

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

Download in other formats:

Original Format