Context Navigation

physiq_mod.F @ 2266

Last change on this file since 2266 was 2265, checked in by emillour, 5 years ago
Mars GCM: Save "dtau", the opacity difference between model and target dust scenario in the restartfi.nc file so that we have 1+1=2 when running with dust injection schemes. EM
File size: 140.7 KB

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

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

Download in other formats:

Original Format