Context Navigation

physiq_mod.F @ 2997

Last change on this file since 2997 was 2965, checked in by abierjon, 20 months ago

Mars GCM:
Fixed r2963 which was preventing to compile the model without XIOS
+ changed the computation of variables rhowater_* so that they are real densities (factor 1/rvap missing ; this doesn't affect the previous PEM results as these densities were only compared between each other)
+ added comments and units for ice table variables in physiq_mod.F
+ made Clapeyron coefficient names in physiq_mod.F coherent with how they are defined in the PEM
+ fixed a reference in constants_marspem_mod.F90
+ fixed unit attribute of surface/soil water densities in field_def_physics_mars.xml

File size: 158.0 KB

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

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

Download in other formats:

Original Format