Context Navigation

physiq_mod.F @ 3568

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

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

Download in other formats:

Original Format