Context Navigation

physiq_mod.F @ 2223

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

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

Download in other formats:

Original Format