Context Navigation

physiq_mod.F @ 2156

Last change on this file since 2156 was 2149, checked in by emillour, 5 years ago
Mars GCM: Add F.Lott's non-orographic GW parametrization. Disabled by default for now, activated by setting calllott_nonoro=.true. in callphys.def GG+EM
File size: 132.5 KB

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

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

Download in other formats:

Original Format