Context Navigation

physiq_mod.F @ 2429

Last change on this file since 2429 was 2417, checked in by emillour, 5 years ago

Mars GCM:
Add a new scheme to handle the dust and its radiative impact. Triggered using
a new flag dustscaling_mode=2 (dustscaling_mod=0: no rescaling at all, and
dustscaling_mode=1: full rescaling using tauscaling, GCMv5.3 style). Rescaling
is then only done on the radiative impact (see dust_scaling_mod.F90) of dust.
Moreover the scaling factor "dust_rad_adjust" is evaluated using the target dust
scenario opacity for the next sol and left to evolve linearly until then to not
impose the diurnal evolution of dust.
In practice, main changes or additions in the code are:

renamed flag "tauscaling_mode" as "dustscaling_mode"
moved parameter "t_scenario_sol" to "dust_param_mod"
adapted "compute_dustscaling" routine in "dust_scaling_mod"
added module "dust_rad_adjust_mod"
2D fields "dust_rad_adjust_prev" and "dust_rad_adjust_next" required to compute coefficient "dust_rad_adjust" need to be stored in (re)startfi files

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

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

Download in other formats:

Original Format