Context Navigation

physiq.F @ 1124

Last change on this file since 1124 was 1124, checked in by emillour, 11 years ago

Mars GCM:

Improved 15um cooling routines, with also better handling of errors.

FGG

File size: 93.1 KB

Rev	Line
[234]	1	SUBROUTINE physiq(
[226]	2	$ ngrid,nlayer,nq
	3	$ ,firstcall,lastcall
	4	$ ,pday,ptime,ptimestep
	5	$ ,pplev,pplay,pphi
	6	$ ,pu,pv,pt,pq
	7	$ ,pw
	8	$ ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn
[86]	9	#ifdef MESOSCALE
[226]	10	#include "meso_inc/meso_inc_invar.F"
[86]	11	#endif
[226]	12	$ )
[42]	13
[1036]	14	use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2,
	15	& igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
	16	& igcm_ccn_mass, igcm_ccn_number,
	17	& igcm_dust_mass, igcm_dust_number, igcm_h2o2,
	18	& nuice_ref, rho_ice, rho_dust, ref_r0
[1047]	19	use comsoil_h, only: inertiedat, ! soil thermal inertia
	20	& nsoilmx ! number of subsurface layers
	21	use eofdump_mod, only: eofdump
	22	use comgeomfi_h, only: long, lati, area
	23	use comdiurn_h, only: sinlon, coslon, sinlat, coslat
	24	use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
	25	& zthe, z0, albedo_h2o_ice,
	26	& frost_albedo_threshold
	27	use comsaison_h, only: dist_sol, declin, mu0, fract
	28	use slope_mod, only: theta_sl, psi_sl
	29	use conc_mod, only: rnew, cpnew, mmean
[1112]	30
	31	#ifdef MESOSCALE
	32	use comsoil_h, only: mlayer,layer
	33	use surfdat_h, only: z0_default
	34	#endif
	35
[42]	36	IMPLICIT NONE
	37	c=======================================================================
	38	c
	39	c subject:
	40	c --------
	41	c
	42	c Organisation of the physical parametrisations of the LMD
	43	c martian atmospheric general circulation model.
	44	c
	45	c The GCM can be run without or with tracer transport
	46	c depending on the value of Logical "tracer" in file "callphys.def"
	47	c Tracers may be water vapor, ice OR chemical species OR dust particles
	48	c
	49	c SEE comments in initracer.F about numbering of tracer species...
	50	c
	51	c It includes:
	52	c
	53	c 1. Initialization:
	54	c 1.1 First call initializations
	55	c 1.2 Initialization for every call to physiq
	56	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	57	c 2. Compute radiative transfer tendencies
	58	c (longwave and shortwave) for CO2 and aerosols.
	59	c 3. Gravity wave and subgrid scale topography drag :
	60	c 4. Vertical diffusion (turbulent mixing):
	61	c 5. Convective adjustment
	62	c 6. Condensation and sublimation of carbon dioxide.
	63	c 7. TRACERS :
	64	c 7a. water and water ice
	65	c 7b. call for photochemistry when tracers are chemical species
	66	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
	67	c 7d. updates (CO2 pressure variations, surface budget)
	68	c 8. Contribution to tendencies due to thermosphere
	69	c 9. Surface and sub-surface temperature calculations
	70	c 10. Write outputs :
	71	c - "startfi", "histfi" (if it's time)
	72	c - Saving statistics (if "callstats = .true.")
	73	c - Dumping eof (if "calleofdump = .true.")
	74	c - Output any needed variables in "diagfi"
	75	c 11. Diagnostic: mass conservation of tracers
	76	c
	77	c author:
	78	c -------
	79	c Frederic Hourdin 15/10/93
	80	c Francois Forget 1994
	81	c Christophe Hourdin 02/1997
	82	c Subroutine completly rewritten by F.Forget (01/2000)
	83	c Introduction of the photochemical module: S. Lebonnois (11/2002)
	84	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
	85	c Water ice clouds: Franck Montmessin (update 06/2003)
	86	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	87	c Nb: See callradite.F for more information.
[234]	88	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
[414]	89	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
[42]	90	c
	91	c arguments:
	92	c ----------
	93	c
	94	c input:
	95	c ------
	96	c ecri period (in dynamical timestep) to write output
	97	c ngrid Size of the horizontal grid.
	98	c All internal loops are performed on that grid.
	99	c nlayer Number of vertical layers.
	100	c nq Number of advected fields
	101	c firstcall True at the first call
	102	c lastcall True at the last call
	103	c pday Number of days counted from the North. Spring
	104	c equinoxe.
	105	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	106	c ptimestep timestep (s)
	107	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	108	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	109	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	110	c pu(ngrid,nlayer) u component of the wind (ms-1)
	111	c pv(ngrid,nlayer) v component of the wind (ms-1)
	112	c pt(ngrid,nlayer) Temperature (K)
	113	c pq(ngrid,nlayer,nq) Advected fields
[330]	114	c pudyn(ngrid,nlayer) \|
	115	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
	116	c ptdyn(ngrid,nlayer) \| corresponding variables
	117	c pqdyn(ngrid,nlayer,nq) \|
[42]	118	c pw(ngrid,?) vertical velocity
	119	c
	120	c output:
	121	c -------
	122	c
[1047]	123	c pdu(ngrid,nlayer) \|
	124	c pdv(ngrid,nlayer) \| Temporal derivative of the corresponding
	125	c pdt(ngrid,nlayer) \| variables due to physical processes.
	126	c pdq(ngrid,nlayer,nq) \|
	127	c pdpsrf(ngrid) \|
[42]	128	c tracerdyn call tracer in dynamical part of GCM ?
	129
	130	c
	131	c=======================================================================
	132	c
	133	c 0. Declarations :
	134	c ------------------
	135
	136	#include "dimensions.h"
	137	#include "dimphys.h"
[1047]	138	!#include "comgeomfi.h"
	139	!#include "surfdat.h"
	140	!#include "comsoil.h"
	141	!#include "comdiurn.h"
[42]	142	#include "callkeys.h"
	143	#include "comcstfi.h"
	144	#include "planete.h"
[1047]	145	!#include "comsaison.h"
[42]	146	#include "control.h"
[1047]	147	!#include "dimradmars.h"
	148	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
	149	#include"scatterers.h"
[42]	150	#include "comg1d.h"
[1036]	151	!#include "tracer.h"
[42]	152	#include "nlteparams.h"
[883]	153	#include "comvert.h"
[42]	154
	155	#include "chimiedata.h"
	156	#include "param.h"
[635]	157	#include "param_v4.h"
[1047]	158	!#include "conc.h"
[42]	159
	160	#include "netcdf.inc"
	161
[1047]	162	!#include "slope.h"
[234]	163
[86]	164	#ifdef MESOSCALE
[42]	165	#include "wrf_output_2d.h"
	166	#include "wrf_output_3d.h"
[1036]	167	!#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
[226]	168	#include "meso_inc/meso_inc_var.F"
[86]	169	#endif
[42]	170
	171	c Arguments :
	172	c -----------
	173
	174	c inputs:
	175	c -------
	176	INTEGER ngrid,nlayer,nq
	177	REAL ptimestep
[1047]	178	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
	179	REAL pphi(ngrid,nlayer)
	180	REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
	181	REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
	182	REAL pw(ngrid,nlayer) !Mars pvervel transmit par dyn3d
	183	REAL zh(ngrid,nlayer) ! potential temperature (K)
[42]	184	LOGICAL firstcall,lastcall
[226]	185
[42]	186	REAL pday
[999]	187	REAL ptime
[42]	188	logical tracerdyn
	189
	190	c outputs:
	191	c --------
	192	c physical tendencies
[1047]	193	REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
	194	REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
	195	REAL pdpsrf(ngrid) ! surface pressure tendency
[42]	196
	197
	198	c Local saved variables:
	199	c ----------------------
	200	c aerosol (dust or ice) extinction optical depth at reference wavelength
[1047]	201	c "longrefvis" set in dimradmars_mod , for one of the "naerkind" kind of
[42]	202	c aerosol optical properties :
[1047]	203	REAL,SAVE,ALLOCATABLE :: aerosol(:,:,:)
[42]	204
[1047]	205	INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
	206	INTEGER,SAVE :: icount ! counter of calls to physiq during the run.
	207	REAL,SAVE,ALLOCATABLE :: tsurf(:) ! Surface temperature (K)
	208	REAL,SAVE,ALLOCATABLE :: tsoil(:,:) ! sub-surface temperatures (K)
	209	REAL,SAVE,ALLOCATABLE :: co2ice(:) ! co2 ice surface layer (kg.m-2)
	210	REAL,SAVE,ALLOCATABLE :: albedo(:,:) ! Surface albedo in each solar band
	211	REAL,SAVE,ALLOCATABLE :: emis(:) ! Thermal IR surface emissivity
	212	REAL,SAVE,ALLOCATABLE :: dtrad(:,:) ! Net atm. radiative heating rate (K.s-1)
	213	REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:) ! rad. flux from sky absorbed by surface (W.m-2)
	214	REAL,SAVE,ALLOCATABLE :: fluxrad(:) ! Net radiative surface flux (W.m-2)
	215	REAL,SAVE,ALLOCATABLE :: capcal(:) ! surface heat capacity (J m-2 K-1)
	216	REAL,SAVE,ALLOCATABLE :: fluxgrd(:) ! surface conduction flux (W.m-2)
[1036]	217	REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
[1047]	218	REAL,SAVE,ALLOCATABLE :: q2(:,:) ! Turbulent Kinetic Energy
[283]	219
[42]	220	c Variables used by the water ice microphysical scheme:
[1047]	221	REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m)
	222	REAL nuice(ngrid,nlayer) ! Estimated effective variance
[42]	223	! of the size distribution
[1047]	224	real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
	225	real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3)
	226	REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
	227	REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry)
	228	REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
	229	! thermal inertia (J.s-1/2.m-2.K-1)
	230	! (used only when tifeedback=.true.)
[42]	231
[234]	232	c Variables used by the slope model
	233	REAL sl_ls, sl_lct, sl_lat
	234	REAL sl_tau, sl_alb, sl_the, sl_psi
	235	REAL sl_fl0, sl_flu
	236	REAL sl_ra, sl_di0
	237	REAL sky
	238
[1047]	239	REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
[42]	240
	241	c Local variables :
	242	c -----------------
	243
	244	REAL CBRT
	245	EXTERNAL CBRT
	246
[1036]	247	! CHARACTER*80 fichier
	248	INTEGER l,ig,ierr,igout,iq,tapphys
[42]	249
[1047]	250	REAL fluxsurf_lw(ngrid) !incident LW (IR) surface flux (W.m-2)
	251	REAL fluxsurf_sw(ngrid,2) !incident SW (solar) surface flux (W.m-2)
	252	REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2)
	253	REAL fluxtop_sw(ngrid,2) !Outgoing SW (solar) flux to space (W.m-2)
	254	REAL tauref(ngrid) ! Reference column optical depth at odpref
[627]	255	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
[1047]	256	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
[42]	257	REAL zls ! solar longitude (rad)
	258	REAL zday ! date (time since Ls=0, in martian days)
[1047]	259	REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers
	260	REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
[1036]	261	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
[42]	262
	263	c Tendancies due to various processes:
[1047]	264	REAL dqsurf(ngrid,nq)
	265	REAL zdtlw(ngrid,nlayer) ! (K/s)
	266	REAL zdtsw(ngrid,nlayer) ! (K/s)
	267	! REAL cldtlw(ngrid,nlayer) ! (K/s) LW heating rate for clear area
	268	! REAL cldtsw(ngrid,nlayer) ! (K/s) SW heating rate for clear area
	269	REAL zdtnirco2(ngrid,nlayer) ! (K/s)
	270	REAL zdtnlte(ngrid,nlayer) ! (K/s)
	271	REAL zdtsurf(ngrid) ! (K/s)
	272	REAL zdtcloud(ngrid,nlayer)
	273	REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
	274	REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid) ! (K/s)
	275	REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
	276	REAL zdhadj(ngrid,nlayer) ! (K/s)
	277	REAL zdtgw(ngrid,nlayer) ! (K/s)
	278	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
	279	REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
	280	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
[42]	281
[1047]	282	REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
	283	REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
	284	REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
	285	REAL zdqadj(ngrid,nlayer,nq)
	286	REAL zdqc(ngrid,nlayer,nq)
	287	REAL zdqcloud(ngrid,nlayer,nq)
	288	REAL zdqscloud(ngrid,nq)
	289	REAL zdqchim(ngrid,nlayer,nq)
	290	REAL zdqschim(ngrid,nq)
[42]	291
[1047]	292	REAL zdteuv(ngrid,nlayer) ! (K/s)
	293	REAL zdtconduc(ngrid,nlayer) ! (K/s)
	294	REAL zdumolvis(ngrid,nlayer)
	295	REAL zdvmolvis(ngrid,nlayer)
	296	real zdqmoldiff(ngrid,nlayer,nq)
[42]	297
	298	c Local variable for local intermediate calcul:
[1047]	299	REAL zflubid(ngrid)
	300	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
	301	REAL zdum1(ngrid,nlayer)
	302	REAL zdum2(ngrid,nlayer)
[42]	303	REAL ztim1,ztim2,ztim3, z1,z2
	304	REAL ztime_fin
[1047]	305	REAL zdh(ngrid,nlayer)
[42]	306	INTEGER length
	307	PARAMETER (length=100)
	308
	309	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	310	c -----------------------------------------------------------------------------
[1047]	311	REAL ps(ngrid), zt(ngrid,nlayer)
	312	REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
	313	REAL zq(ngrid,nlayer,nq)
	314	REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
[42]	315	character*2 str2
[1036]	316	! character*5 str5
[1047]	317	real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
	318	REAL,SAVE,ALLOCATABLE :: tauscaling(:) ! Convertion factor for qdust and Ndust
	319	real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
[42]	320	integer igmin, lmin
	321	logical tdiag
	322
[1047]	323	real co2col(ngrid) ! CO2 column
[883]	324	! pplev and pplay are dynamical inputs and must not be modified in the physics.
	325	! instead, use zplay and zplev :
[1047]	326	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
[1036]	327	! REAL zstress(ngrid),cd
[1047]	328	real tmean, zlocal(nlayer)
	329	real rho(ngrid,nlayer) ! density
	330	real vmr(ngrid,nlayer) ! volume mixing ratio
	331	real rhopart(ngrid,nlayer) ! number density of a given species
	332	real colden(ngrid,nq) ! vertical column of tracers
	333	REAL mtot(ngrid) ! Total mass of water vapor (kg/m2)
	334	REAL icetot(ngrid) ! Total mass of water ice (kg/m2)
	335	REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2)
	336	REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2)
	337	REAL rave(ngrid) ! Mean water ice effective radius (m)
	338	REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1
	339	REAL tauTES(ngrid) ! column optical depth at 825 cm-1
[42]	340	REAL Qabsice ! Water ice absorption coefficient
[1047]	341	REAL taucloudtes(ngrid) ! Cloud opacity at infrared
[520]	342	! reference wavelength using
	343	! Qabs instead of Qext
	344	! (direct comparison with TES)
[719]	345
[1047]	346	REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
	347	REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
	348	REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
	349	REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
	350	REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg)
	351	REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg)
[42]	352
[411]	353	c Test 1d/3d scavenging
[1047]	354	real h2otot(ngrid)
	355	REAL satu(ngrid,nlayer) ! satu ratio for output
	356	REAL zqsat(ngrid,nlayer) ! saturation
[42]	357
[1047]	358	REAL,SAVE :: time_phys
[42]	359
[414]	360	! Added for new NLTE scheme
	361
[1047]	362	real co2vmr_gcm(ngrid,nlayer)
	363	real n2vmr_gcm(ngrid,nlayer)
	364	real ovmr_gcm(ngrid,nlayer)
	365	real covmr_gcm(ngrid,nlayer)
[1124]	366	integer ierr_nlte
	367	real*8 varerr
[414]	368
[267]	369	c Variables for PBL
[1047]	370	REAL zz1(ngrid)
	371	REAL lmax_th_out(ngrid),zmax_th(ngrid)
	372	REAL,SAVE,ALLOCATABLE :: wstar(:)
	373	REAL,SAVE,ALLOCATABLE :: hfmax_th(:)
	374	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
	375	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
	376	INTEGER lmax_th(ngrid),dimout,n_out,n
[566]	377	CHARACTER(50) zstring
[1047]	378	REAL dtke_th(ngrid,nlayer+1)
	379	REAL zcdv(ngrid), zcdh(ngrid)
[636]	380	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
[566]	381	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
[1047]	382	REAL u_out1(ngrid)
	383	REAL T_out1(ngrid)
[566]	384	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
[1047]	385	REAL ustar(ngrid),tstar(ngrid) ! friction velocity and friction potential temp
	386	REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
	387	! REAL zu2(ngrid)
	388	REAL sensibFlux(ngrid)
[790]	389
[42]	390	c=======================================================================
	391
	392	c 1. Initialisation:
	393	c -----------------
	394
	395	c 1.1 Initialisation only at first call
	396	c ---------------------------------------
	397	IF (firstcall) THEN
	398
[1036]	399	! allocate local (saved) arrays:
	400	allocate(qsurf(ngrid,nq))
[1047]	401	allocate(tsoil(ngrid,nsoilmx))
	402	allocate(tsurf(ngrid))
	403	allocate(aerosol(ngrid,nlayer,naerkind))
	404	allocate(co2ice(ngrid))
	405	allocate(albedo(ngrid,2))
	406	allocate(emis(ngrid))
	407	allocate(dtrad(ngrid,nlayer))
	408	allocate(fluxrad_sky(ngrid))
	409	allocate(fluxrad(ngrid))
	410	allocate(capcal(ngrid))
	411	allocate(fluxgrd(ngrid))
	412	allocate(q2(ngrid,nlayer+1))
	413	allocate(tauscaling(ngrid))
	414	allocate(wstar(ngrid))
	415	allocate(hfmax_th(ngrid))
[1036]	416
[42]	417	c variables set to 0
	418	c ~~~~~~~~~~~~~~~~~~
[286]	419	aerosol(:,:,:)=0
	420	dtrad(:,:)=0
[674]	421
	422	#ifndef MESOSCALE
[286]	423	fluxrad(:)=0
[528]	424	wstar(:)=0.
[674]	425	#else
	426	#include "meso_inc/meso_inc_ini_restart.F"
	427	#endif
[268]	428
[42]	429	c read startfi
	430	c ~~~~~~~~~~~~
[226]	431	#ifndef MESOSCALE
	432	! Read netcdf initial physical parameters.
	433	CALL phyetat0 ("startfi.nc",0,0,
[1047]	434	& nsoilmx,ngrid,nlayer,nq,
[226]	435	& day_ini,time_phys,
	436	& tsurf,tsoil,emis,q2,qsurf,co2ice)
	437	#else
	438	#include "meso_inc/meso_inc_ini.F"
	439	#endif
[185]	440
[42]	441	if (pday.ne.day_ini) then
	442	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	443	& "physics and dynamics"
	444	write(,) "dynamics day: ",pday
	445	write(,) "physics day: ",day_ini
	446	stop
	447	endif
	448
	449	write (,) 'In physiq day_ini =', day_ini
	450
[286]	451	c initialize tracers
	452	c ~~~~~~~~~~~~~~~~~~
	453	tracerdyn=tracer
	454	IF (tracer) THEN
[1036]	455	CALL initracer(ngrid,nq,qsurf,co2ice)
[286]	456	ENDIF ! end tracer
	457
[42]	458	c Initialize albedo and orbital calculation
	459	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	460	CALL surfini(ngrid,co2ice,qsurf,albedo)
	461	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	462
	463	c initialize soil
	464	c ~~~~~~~~~~~~~~~
	465	IF (callsoil) THEN
[833]	466	c Thermal inertia feedback:
	467	IF (tifeedback) THEN
	468	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
	469	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
	470	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	471	ELSE
	472	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	473	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	474	ENDIF ! of IF (tifeedback)
[42]	475	ELSE
	476	PRINT*,
	477	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	478	DO ig=1,ngrid
	479	capcal(ig)=1.e5
	480	fluxgrd(ig)=0.
	481	ENDDO
	482	ENDIF
	483	icount=1
	484
[226]	485	#ifndef MESOSCALE
	486	c Initialize thermospheric parameters
	487	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[42]	488
[705]	489	if (callthermos) then
	490	if(solvarmod.eq.0) call param_read
	491	if(solvarmod.eq.1) call param_read_e107
	492	endif
[226]	493	#endif
[42]	494	c Initialize R and Cp as constant
	495
	496	if (.not.callthermos .and. .not.photochem) then
[1047]	497	do l=1,nlayer
	498	do ig=1,ngrid
[42]	499	rnew(ig,l)=r
	500	cpnew(ig,l)=cpp
	501	mmean(ig,l)=mugaz
	502	enddo
	503	enddo
	504	endif
	505
[757]	506	if(callnlte.and.nltemodel.eq.2) call nlte_setup
[414]	507	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
[1013]	508	if(thermochem) call chemthermos_readini
[414]	509
[1047]	510	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
[283]	511	write(,)"physiq: water_param Surface water ice albedo:",
	512	. albedo_h2o_ice
[42]	513	ENDIF
[900]	514
	515	#ifndef MESOSCALE
[1047]	516	if (callslope) call getslopes(ngrid,phisfi)
[999]	517
[1047]	518	call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq,
[999]	519	. ptimestep,pday,time_phys,area,
	520	. albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
	521
[900]	522	#endif
[42]	523
	524	ENDIF ! (end of "if firstcall")
	525
[414]	526
[42]	527	c ---------------------------------------------------
	528	c 1.2 Initializations done at every physical timestep:
	529	c ---------------------------------------------------
	530	c
	531
	532	c Initialize various variables
	533	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[286]	534	pdv(:,:)=0
	535	pdu(:,:)=0
	536	pdt(:,:)=0
	537	pdq(:,:,:)=0
	538	pdpsrf(:)=0
	539	zflubid(:)=0
	540	zdtsurf(:)=0
	541	dqsurf(:,:)=0
[42]	542	igout=ngrid/2+1
	543
	544
	545	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	546
	547	c Compute Solar Longitude (Ls) :
	548	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	549	if (season) then
	550	call solarlong(zday,zls)
	551	else
	552	call solarlong(float(day_ini),zls)
	553	end if
	554
[883]	555	c Initialize pressure levels
	556	c ~~~~~~~~~~~~~~~~~~
	557	zplev(:,:) = pplev(:,:)
	558	zplay(:,:) = pplay(:,:)
	559	ps(:) = pplev(:,1)
	560
	561
[42]	562	c Compute geopotential at interlayers
	563	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	564	c ponderation des altitudes au niveau des couches en dp/p
	565
	566	DO l=1,nlayer
	567	DO ig=1,ngrid
	568	zzlay(ig,l)=pphi(ig,l)/g
	569	ENDDO
	570	ENDDO
	571	DO ig=1,ngrid
	572	zzlev(ig,1)=0.
	573	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	574	ENDDO
	575	DO l=2,nlayer
	576	DO ig=1,ngrid
[883]	577	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
	578	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
[42]	579	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	580	ENDDO
	581	ENDDO
	582
	583
	584	! Potential temperature calculation not the same in physiq and dynamic
	585
	586	c Compute potential temperature
	587	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	588	DO l=1,nlayer
	589	DO ig=1,ngrid
[883]	590	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
[42]	591	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	592	ENDDO
	593	ENDDO
	594
[226]	595	#ifndef MESOSCALE
	596	c-----------------------------------------------------------------------
	597	c 1.2.5 Compute mean mass, cp, and R
	598	c --------------------------------
[42]	599
[226]	600	if(photochem.or.callthermos) then
[1047]	601	call concentrations(ngrid,nlayer,nq,
	602	& zplay,pt,pdt,pq,pdq,ptimestep)
[226]	603	endif
	604	#endif
[42]	605	c-----------------------------------------------------------------------
	606	c 2. Compute radiative tendencies :
	607	c------------------------------------
	608
	609
	610	IF (callrad) THEN
	611	IF( MOD(icount-1,iradia).EQ.0) THEN
	612
	613	c Local Solar zenith angle
	614	c ~~~~~~~~~~~~~~~~~~~~~~~~
	615	CALL orbite(zls,dist_sol,declin)
	616
	617	IF(diurnal) THEN
	618	ztim1=SIN(declin)
	619	ztim2=COS(declin)COS(2.pi*(zday-.5))
	620	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	621
	622	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	623	s ztim1,ztim2,ztim3, mu0,fract)
	624
	625	ELSE
	626	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
	627	ENDIF
	628
	629	c NLTE cooling from CO2 emission
	630	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[414]	631	IF(callnlte) then
	632	if(nltemodel.eq.0.or.nltemodel.eq.1) then
[883]	633	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
[414]	634	else if(nltemodel.eq.2) then
[705]	635	co2vmr_gcm(1:ngrid,1:nlayer)=
	636	& pq(1:ngrid,1:nlayer,igcm_co2)*
	637	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
	638	n2vmr_gcm(1:ngrid,1:nlayer)=
	639	& pq(1:ngrid,1:nlayer,igcm_n2)*
	640	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
	641	covmr_gcm(1:ngrid,1:nlayer)=
	642	& pq(1:ngrid,1:nlayer,igcm_co)*
	643	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
	644	ovmr_gcm(1:ngrid,1:nlayer)=
	645	& pq(1:ngrid,1:nlayer,igcm_o)*
	646	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
[414]	647
[883]	648	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
[414]	649	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
[1124]	650	$ ovmr_gcm, zdtnlte,ierr_nlte,varerr )
	651	if(ierr_nlte.gt.0) then
	652	write(,)
	653	$ 'WARNING: nlte_tcool output with error message',
	654	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
	655	write(,)'I will continue anyway'
	656	endif
[42]	657
[705]	658	zdtnlte(1:ngrid,1:nlayer)=
	659	& zdtnlte(1:ngrid,1:nlayer)/86400.
[414]	660	endif
[528]	661	else
	662	zdtnlte(:,:)=0.
[414]	663	endif
[42]	664
	665	c Find number of layers for LTE radiation calculations
	666	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
[883]	667	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
[42]	668
	669	c Note: Dustopacity.F has been transferred to callradite.F
	670
	671	c Call main radiative transfer scheme
	672	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	673	c Transfer through CO2 (except NIR CO2 absorption)
	674	c and aerosols (dust and water ice)
	675
	676	c Radiative transfer
	677	c ------------------
[833]	678
[42]	679	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
[883]	680	$ emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
[42]	681	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
[520]	682	$ tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
	683	$ nuice,co2ice)
[42]	684
[234]	685	c Outputs for basic check (middle of domain)
	686	c ------------------------------------------
[627]	687	write(*,'("Ls =",f11.6," check lat =",f10.6,
	688	& " lon =",f11.6)')
	689	& zls180./pi,lati(igout)180/pi,long(igout)*180/pi
	690	write(*,'(" tauref(",f4.0," Pa) =",f9.6,
	691	& " tau(",f4.0," Pa) =",f9.6)')
	692	& odpref,tauref(igout),
[883]	693	& odpref,tau(igout,1)*odpref/zplev(igout,1)
[234]	694	c ---------------------------------------------------------
	695	c Call slope parameterization for direct and scattered flux
	696	c ---------------------------------------------------------
	697	IF(callslope) THEN
	698	print *, 'Slope scheme is on and computing...'
	699	DO ig=1,ngrid
	700	sl_the = theta_sl(ig)
	701	IF (sl_the .ne. 0.) THEN
	702	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
	703	DO l=1,2
	704	sl_lct = ptime24. + 180.long(ig)/pi/15.
	705	sl_ra = pi*(1.0-sl_lct/12.)
	706	sl_lat = 180.*lati(ig)/pi
[577]	707	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
[234]	708	sl_alb = albedo(ig,l)
	709	sl_psi = psi_sl(ig)
	710	sl_fl0 = fluxsurf_sw(ig,l)
	711	sl_di0 = 0.
	712	if (mu0(ig) .gt. 0.) then
	713	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
	714	sl_di0 = sl_di0*1370./dist_sol/dist_sol
	715	sl_di0 = sl_di0/ztim1
	716	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
	717	endif
	718	! you never know (roundup concern...)
	719	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
	720	!!!!!!!!!!!!!!!!!!!!!!!!!!
	721	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
	722	& sl_tau, sl_alb, sl_the, sl_psi,
	723	& sl_di0, sl_fl0, sl_flu )
	724	!!!!!!!!!!!!!!!!!!!!!!!!!!
	725	fluxsurf_sw(ig,l) = sl_flu
	726	ENDDO
	727	!!! compute correction on IR flux as well
	728	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
	729	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
	730	ENDIF
	731	ENDDO
	732	ENDIF
	733
[42]	734	c CO2 near infrared absorption
	735	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[286]	736	zdtnirco2(:,:)=0
[42]	737	if (callnirco2) then
[883]	738	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
[42]	739	. mu0,fract,declin, zdtnirco2)
	740	endif
	741
	742	c Radiative flux from the sky absorbed by the surface (W.m-2)
	743	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	744	DO ig=1,ngrid
	745	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	746	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
	747	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
	748	ENDDO
	749
	750
	751	c Net atmospheric radiative heating rate (K.s-1)
	752	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	753	IF(callnlte) THEN
[883]	754	CALL blendrad(ngrid, nlayer, zplay,
[42]	755	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
	756	ELSE
	757	DO l=1,nlayer
	758	DO ig=1,ngrid
	759	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	760	& +zdtnirco2(ig,l)
	761	ENDDO
	762	ENDDO
	763	ENDIF
	764
	765	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	766
	767	c Transformation of the radiative tendencies:
	768	c -------------------------------------------
	769
	770	c Net radiative surface flux (W.m-2)
	771	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	772	c
	773	DO ig=1,ngrid
	774	zplanck(ig)=tsurf(ig)*tsurf(ig)
	775	zplanck(ig)=emis(ig)*
	776	$ stephanzplanck(ig)zplanck(ig)
	777	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
[234]	778	IF(callslope) THEN
	779	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
	780	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
	781	ENDIF
[42]	782	ENDDO
	783
	784	DO l=1,nlayer
	785	DO ig=1,ngrid
	786	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	787	ENDDO
	788	ENDDO
	789
	790	ENDIF ! of IF (callrad)
	791
[226]	792	c-----------------------------------------------------------------------
	793	c 3. Gravity wave and subgrid scale topography drag :
	794	c -------------------------------------------------
[42]	795
[226]	796
	797	IF(calllott)THEN
	798
	799	CALL calldrag_noro(ngrid,nlayer,ptimestep,
[883]	800	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
[226]	801
	802	DO l=1,nlayer
	803	DO ig=1,ngrid
	804	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
	805	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
	806	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
	807	ENDDO
	808	ENDDO
	809	ENDIF
[234]	810
[42]	811	c-----------------------------------------------------------------------
	812	c 4. Vertical diffusion (turbulent mixing):
	813	c -----------------------------------------
[226]	814
[42]	815	IF (calldifv) THEN
	816
	817	DO ig=1,ngrid
	818	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	819	ENDDO
	820
[286]	821	zdum1(:,:)=0
	822	zdum2(:,:)=0
[42]	823	do l=1,nlayer
	824	do ig=1,ngrid
	825	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	826	enddo
	827	enddo
[226]	828
[288]	829
	830	#ifdef MESOSCALE
[289]	831	IF (.not.flag_LES) THEN
[288]	832	#endif
	833	c ----------------------
[284]	834	c Treatment of a special case : using new surface layer (Richardson based)
	835	c without using the thermals in gcm and mesoscale can yield problems in
	836	c weakly unstable situations when winds are near to 0. For those cases, we add
	837	c a unit subgrid gustiness. Remember that thermals should be used we using the
	838	c Richardson based surface layer model.
[288]	839	IF ( .not.calltherm .and. callrichsl ) THEN
[1047]	840	DO ig=1, ngrid
[284]	841	IF (zh(ig,1) .lt. tsurf(ig)) THEN
[528]	842	wstar(ig)=1.
	843	hfmax_th(ig)=0.2
	844	ELSE
	845	wstar(ig)=0.
	846	hfmax_th(ig)=0.
	847	ENDIF
[284]	848	ENDDO
	849	ENDIF
[288]	850	c ----------------------
	851	#ifdef MESOSCALE
[284]	852	ENDIF
	853	#endif
	854
[544]	855	IF (tke_heat_flux .ne. 0.) THEN
	856	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
[883]	857	& (-alog(zplay(:,1)/zplev(:,1)))
[544]	858	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
	859	ENDIF
[288]	860
[42]	861	c Calling vdif (Martian version WITH CO2 condensation)
	862	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
	863	$ ptimestep,capcal,lwrite,
[883]	864	$ zplay,zplev,zzlay,zzlev,z0,
[42]	865	$ pu,pv,zh,pq,tsurf,emis,qsurf,
	866	$ zdum1,zdum2,zdh,pdq,zflubid,
	867	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
[660]	868	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux
[529]	869	#ifdef MESOSCALE
	870	& ,flag_LES
	871	#endif
	872	& )
[42]	873
[529]	874
[226]	875	#ifdef MESOSCALE
	876	#include "meso_inc/meso_inc_les.F"
	877	#endif
[42]	878	DO l=1,nlayer
	879	DO ig=1,ngrid
	880	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	881	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	882	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	883
	884	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	885	ENDDO
	886	ENDDO
	887
[226]	888	DO ig=1,ngrid
	889	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	890	ENDDO
[42]	891
[790]	892	if (tracer) then
	893	#ifdef MESOSCALE
	894	IF (.not.flag_LES) THEN
	895	#endif
[42]	896	DO iq=1, nq
	897	DO l=1,nlayer
	898	DO ig=1,ngrid
	899	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	900	ENDDO
	901	ENDDO
	902	ENDDO
	903	DO iq=1, nq
	904	DO ig=1,ngrid
	905	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	906	ENDDO
	907	ENDDO
[790]	908	#ifdef MESOSCALE
	909	ENDIF
	910	#endif
[42]	911	end if ! of if (tracer)
	912
	913	ELSE
	914	DO ig=1,ngrid
	915	zdtsurf(ig)=zdtsurf(ig)+
	916	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	917	ENDDO
[226]	918	#ifdef MESOSCALE
[42]	919	IF (flag_LES) THEN
	920	write(,) 'LES mode !'
	921	write(,) 'Please set calldifv to T in callphys.def'
	922	STOP
	923	ENDIF
[226]	924	#endif
[42]	925	ENDIF ! of IF (calldifv)
	926
[162]	927	c-----------------------------------------------------------------------
[566]	928	c 5. Thermals :
[162]	929	c -----------------------------
[566]	930
[162]	931	if(calltherm) then
	932
[1032]	933	call calltherm_interface(ngrid,nlayer,nq,
	934	$ tracer,igcm_co2,
[652]	935	$ zzlev,zzlay,
[162]	936	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
[883]	937	$ zplay,zplev,pphi,zpopsk,
[185]	938	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
[660]	939	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
[633]	940
[162]	941	DO l=1,nlayer
	942	DO ig=1,ngrid
	943	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
	944	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
	945	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
	946	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
	947	ENDDO
	948	ENDDO
	949
	950	DO ig=1,ngrid
	951	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
	952	ENDDO
	953
	954	if (tracer) then
	955	DO iq=1,nq
	956	DO l=1,nlayer
	957	DO ig=1,ngrid
	958	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
	959	ENDDO
	960	ENDDO
	961	ENDDO
	962	endif
[42]	963
[277]	964	lmax_th_out(:)=real(lmax_th(:))
	965
[162]	966	else !of if calltherm
	967	lmax_th(:)=0
[528]	968	wstar(:)=0.
	969	hfmax_th(:)=0.
[277]	970	lmax_th_out(:)=0.
[162]	971	end if
	972
[42]	973	c-----------------------------------------------------------------------
	974	c 5. Dry convective adjustment:
	975	c -----------------------------
	976
	977	IF(calladj) THEN
	978
	979	DO l=1,nlayer
	980	DO ig=1,ngrid
	981	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	982	ENDDO
	983	ENDDO
[286]	984	zduadj(:,:)=0
	985	zdvadj(:,:)=0
	986	zdhadj(:,:)=0
[42]	987
	988	CALL convadj(ngrid,nlayer,nq,ptimestep,
[883]	989	$ zplay,zplev,zpopsk,lmax_th,
[42]	990	$ pu,pv,zh,pq,
	991	$ pdu,pdv,zdh,pdq,
	992	$ zduadj,zdvadj,zdhadj,
	993	$ zdqadj)
	994
[162]	995
[42]	996	DO l=1,nlayer
	997	DO ig=1,ngrid
	998	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	999	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	1000	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	1001
	1002	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	1003	ENDDO
	1004	ENDDO
	1005
	1006	if(tracer) then
	1007	DO iq=1, nq
	1008	DO l=1,nlayer
	1009	DO ig=1,ngrid
	1010	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	1011	ENDDO
	1012	ENDDO
	1013	ENDDO
	1014	end if
	1015	ENDIF ! of IF(calladj)
	1016
	1017
[332]	1018
[42]	1019	c-----------------------------------------------------------------------
[1114]	1020	c 6. Specific parameterizations for tracers
[42]	1021	c: -----------------------------------------
	1022
	1023	if (tracer) then
	1024
[1114]	1025	c 6a. Water and ice
[42]	1026	c ---------------
	1027
	1028	c ---------------------------------------
	1029	c Water ice condensation in the atmosphere
	1030	c ----------------------------------------
	1031	IF (water) THEN
	1032
	1033	call watercloud(ngrid,nlayer,ptimestep,
[883]	1034	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
[626]	1035	& pq,pdq,zdqcloud,zdtcloud,
[358]	1036	& nq,tau,tauscaling,rdust,rice,nuice,
	1037	& rsedcloud,rhocloud)
[833]	1038
[633]	1039	c Temperature variation due to latent heat release
[42]	1040	if (activice) then
[633]	1041	pdt(1:ngrid,1:nlayer) =
	1042	& pdt(1:ngrid,1:nlayer) +
	1043	& zdtcloud(1:ngrid,1:nlayer)
[42]	1044	endif
[633]	1045
[42]	1046	! increment water vapour and ice atmospheric tracers tendencies
[706]	1047	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
	1048	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
	1049	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
	1050	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
	1051	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
	1052	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
	1053
	1054	! increment dust and ccn masses and numbers
[883]	1055	! We need to check that we have Nccn & Ndust > 0
	1056	! This is due to single precision rounding problems
[706]	1057	if (microphys) then
	1058	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
	1059	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
	1060	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
	1061	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
	1062	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
	1063	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
[883]	1064	where (pq(:,:,igcm_ccn_mass) +
	1065	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
	1066	pdq(:,:,igcm_ccn_mass) =
	1067	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
	1068	pdq(:,:,igcm_ccn_number) =
	1069	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1070	end where
	1071	where (pq(:,:,igcm_ccn_number) +
	1072	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
	1073	pdq(:,:,igcm_ccn_mass) =
	1074	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
	1075	pdq(:,:,igcm_ccn_number) =
	1076	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1077	end where
[706]	1078	endif
	1079
[883]	1080	if (scavenging) then
	1081	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
	1082	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
	1083	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
	1084	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
	1085	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
	1086	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
	1087	where (pq(:,:,igcm_dust_mass) +
	1088	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
	1089	pdq(:,:,igcm_dust_mass) =
	1090	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
	1091	pdq(:,:,igcm_dust_number) =
	1092	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1093	end where
	1094	where (pq(:,:,igcm_dust_number) +
	1095	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
	1096	pdq(:,:,igcm_dust_mass) =
	1097	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
	1098	pdq(:,:,igcm_dust_number) =
	1099	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1100	end where
	1101	endif ! of if scavenging
[633]	1102
	1103
[42]	1104	END IF ! of IF (water)
	1105
[1114]	1106	c 6b. Aerosol particles
[42]	1107	c -------------------
	1108
	1109	c ----------
	1110	c Dust devil :
	1111	c ----------
	1112	IF(callddevil) then
[883]	1113	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
[42]	1114	& zdqdev,zdqsdev)
	1115
	1116	if (dustbin.ge.1) then
	1117	do iq=1,nq
	1118	DO l=1,nlayer
	1119	DO ig=1,ngrid
	1120	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	1121	ENDDO
	1122	ENDDO
	1123	enddo
	1124	do iq=1,nq
	1125	DO ig=1,ngrid
	1126	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
	1127	ENDDO
	1128	enddo
	1129	endif ! of if (dustbin.ge.1)
	1130
	1131	END IF ! of IF (callddevil)
	1132
	1133	c -------------
	1134	c Sedimentation : acts also on water ice
	1135	c -------------
	1136	IF (sedimentation) THEN
	1137	!call zerophys(ngridnlayernq, zdqsed)
	1138	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	1139	!call zerophys(ngrid*nq, zdqssed)
	1140	zdqssed(1:ngrid,1:nq)=0
	1141
	1142	call callsedim(ngrid,nlayer, ptimestep,
[1005]	1143	& zplev,zzlev, zzlay, pt, pdt, rdust, rice,
[358]	1144	& rsedcloud,rhocloud,
[411]	1145	& pq, pdq, zdqsed, zdqssed,nq,
	1146	& tau,tauscaling)
[833]	1147
[42]	1148	DO iq=1, nq
	1149	DO l=1,nlayer
	1150	DO ig=1,ngrid
	1151	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1152	ENDDO
	1153	ENDDO
	1154	ENDDO
	1155	DO iq=1, nq
	1156	DO ig=1,ngrid
	1157	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1158	ENDDO
	1159	ENDDO
	1160	END IF ! of IF (sedimentation)
[790]	1161
	1162	c Add lifted dust to tendancies after sedimentation in the LES
	1163	#ifdef MESOSCALE
	1164	#include "meso_inc/meso_inc_lift_les.F"
	1165	#endif
	1166
[556]	1167	c
[1114]	1168	c 6c. Chemical species
[556]	1169	c ------------------
[42]	1170
[556]	1171	#ifndef MESOSCALE
	1172	c --------------
	1173	c photochemistry :
	1174	c --------------
	1175	IF (photochem .or. thermochem) then
[411]	1176
[556]	1177	! dust and ice surface area
[883]	1178	call surfacearea(ngrid, nlayer, ptimestep, zplay, zzlay,
[556]	1179	$ pt, pq, pdq, nq,
	1180	$ rdust, rice, tau, tauscaling,
	1181	$ surfdust, surfice)
	1182	! call photochemistry
[1047]	1183	call calchim(ngrid,nlayer,nq,
[1036]	1184	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
[556]	1185	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
	1186	$ zdqcloud,zdqscloud,tauref,co2ice,
	1187	$ pu,pdu,pv,pdv,surfdust,surfice)
	1188
	1189	! increment values of tracers:
	1190	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1191	! tracers is zero anyways
	1192	DO l=1,nlayer
	1193	DO ig=1,ngrid
	1194	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	1195	ENDDO
	1196	ENDDO
	1197	ENDDO ! of DO iq=1,nq
	1198
	1199	! add condensation tendency for H2O2
	1200	if (igcm_h2o2.ne.0) then
	1201	DO l=1,nlayer
	1202	DO ig=1,ngrid
	1203	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	1204	& +zdqcloud(ig,l,igcm_h2o2)
	1205	ENDDO
	1206	ENDDO
	1207	endif
	1208
	1209	! increment surface values of tracers:
	1210	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1211	! tracers is zero anyways
	1212	DO ig=1,ngrid
	1213	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	1214	ENDDO
	1215	ENDDO ! of DO iq=1,nq
	1216
	1217	! add condensation tendency for H2O2
	1218	if (igcm_h2o2.ne.0) then
	1219	DO ig=1,ngrid
	1220	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	1221	& +zdqscloud(ig,igcm_h2o2)
	1222	ENDDO
	1223	endif
	1224
	1225	END IF ! of IF (photochem.or.thermochem)
	1226	#endif
	1227
[1114]	1228	c 6d. Updates
[42]	1229	c ---------
	1230
	1231	DO iq=1, nq
	1232	DO ig=1,ngrid
	1233
	1234	c ---------------------------------
	1235	c Updating tracer budget on surface
	1236	c ---------------------------------
	1237	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1238
	1239	ENDDO ! (ig)
	1240	ENDDO ! (iq)
	1241
	1242	endif ! of if (tracer)
	1243
[226]	1244	#ifndef MESOSCALE
	1245	c-----------------------------------------------------------------------
[1114]	1246	c 7. THERMOSPHERE CALCULATION
[226]	1247	c-----------------------------------------------------------------------
[42]	1248
[226]	1249	if (callthermos) then
[1047]	1250	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
[226]	1251	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	1252	& pt,pq,pu,pv,pdt,pdq,
	1253	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	1254
	1255	DO l=1,nlayer
	1256	DO ig=1,ngrid
	1257	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	1258	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	1259	& +zdteuv(ig,l)
	1260	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	1261	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	1262	DO iq=1, nq
	1263	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	1264	ENDDO
	1265	ENDDO
	1266	ENDDO
	1267
	1268	endif ! of if (callthermos)
	1269	#endif
	1270
[42]	1271	c-----------------------------------------------------------------------
[1114]	1272	c 8. Carbon dioxide condensation-sublimation:
	1273	c (should be the last atmospherical physical process to be computed)
	1274	c -------------------------------------------
	1275
	1276	IF (tituscap) THEN
	1277	!!! get the actual co2 seasonal cap from Titus observations
	1278	CALL geticecover( ngrid, 180.*zls/pi,
	1279	. 180.long/pi, 180.lati/pi, co2ice )
	1280	co2ice = co2ice * 10000.
	1281	ENDIF
	1282
	1283
	1284	pdpsrf(:) = 0
	1285
	1286	IF (callcond) THEN
	1287	CALL newcondens(ngrid,nlayer,nq,ptimestep,
	1288	$ capcal,zplay,zplev,tsurf,pt,
	1289	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	1290	$ co2ice,albedo,emis,
	1291	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
	1292	$ fluxsurf_sw,zls)
	1293
	1294	DO l=1,nlayer
	1295	DO ig=1,ngrid
	1296	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	1297	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	1298	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	1299	ENDDO
	1300	ENDDO
	1301	DO ig=1,ngrid
	1302	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	1303	ENDDO
	1304
	1305	IF (tracer) THEN
	1306	DO iq=1, nq
	1307	DO l=1,nlayer
	1308	DO ig=1,ngrid
	1309	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	1310	ENDDO
	1311	ENDDO
	1312	ENDDO
	1313	ENDIF ! of IF (tracer)
	1314
	1315	#ifndef MESOSCALE
	1316	! update surface pressure
	1317	DO ig=1,ngrid
	1318	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
	1319	ENDDO
	1320
	1321	! update pressure levels
	1322	DO l=1,nlayer
	1323	DO ig=1,ngrid
	1324	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
	1325	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
	1326	ENDDO
	1327	ENDDO
	1328	zplev(:,nlayer+1) = 0.
	1329
	1330	! update layers altitude
	1331	DO l=2,nlayer
	1332	DO ig=1,ngrid
	1333	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
	1334	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
	1335	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	1336	ENDDO
	1337	ENDDO
	1338	#endif
	1339
	1340	ENDIF ! of IF (callcond)
	1341
	1342
	1343	c-----------------------------------------------------------------------
[42]	1344	c 9. Surface and sub-surface soil temperature
	1345	c-----------------------------------------------------------------------
	1346	c
	1347	c
	1348	c 9.1 Increment Surface temperature:
	1349	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1350
	1351	DO ig=1,ngrid
	1352	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1353	ENDDO
	1354
	1355	c Prescribe a cold trap at south pole (except at high obliquity !!)
	1356	c Temperature at the surface is set there to be the temperature
	1357	c corresponding to equilibrium temperature between phases of CO2
	1358
[411]	1359
[1047]	1360	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
[1114]	1361	!#ifndef MESOSCALE
	1362	! if (caps.and.(obliquit.lt.27.)) then => now done in newcondens
[226]	1363	! NB: Updated surface pressure, at grid point 'ngrid', is
[883]	1364	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
[1114]	1365	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	1366	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	1367	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
	1368	! endif
	1369	!#endif
[42]	1370	c -------------------------------------------------------------
[283]	1371	c Change of surface albedo in case of ground frost
[42]	1372	c everywhere except on the north permanent cap and in regions
	1373	c covered by dry ice.
	1374	c ALWAYS PLACE these lines after newcondens !!!
	1375	c -------------------------------------------------------------
	1376	do ig=1,ngrid
	1377	if ((co2ice(ig).eq.0).and.
[283]	1378	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
	1379	albedo(ig,1) = albedo_h2o_ice
	1380	albedo(ig,2) = albedo_h2o_ice
	1381	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
	1382	c write(,) "physiq.F frost :"
	1383	c & ,lati(ig)180./pi, long(ig)180./pi
[42]	1384	endif
	1385	enddo ! of do ig=1,ngrid
[1047]	1386	ENDIF ! of IF (tracer.AND.water.AND.(ngrid.NE.1))
[42]	1387
[528]	1388	c
[42]	1389	c 9.2 Compute soil temperatures and subsurface heat flux:
	1390	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1391	IF (callsoil) THEN
[833]	1392	c Thermal inertia feedback
	1393	IF (tifeedback) THEN
	1394	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
	1395	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
	1396	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1397	ELSE
[42]	1398	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
[833]	1399	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1400	ENDIF
[42]	1401	ENDIF
[411]	1402
[42]	1403
	1404	c-----------------------------------------------------------------------
	1405	c 10. Write output files
	1406	c ----------------------
	1407
[674]	1408	c Save variables for eventual restart in MMM and LES
	1409	#ifdef MESOSCALE
	1410	#include "meso_inc/meso_inc_save_restart.F"
	1411	#endif
	1412
[42]	1413	c -------------------------------
	1414	c Dynamical fields incrementation
	1415	c -------------------------------
	1416	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1417	! temperature, zonal and meridional wind
	1418	DO l=1,nlayer
	1419	DO ig=1,ngrid
	1420	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1421	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1422	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1423	ENDDO
	1424	ENDDO
	1425
	1426	! tracers
	1427	DO iq=1, nq
	1428	DO l=1,nlayer
	1429	DO ig=1,ngrid
	1430	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1431	ENDDO
	1432	ENDDO
	1433	ENDDO
	1434
	1435	! Density
	1436	DO l=1,nlayer
	1437	DO ig=1,ngrid
	1438	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1439	ENDDO
	1440	ENDDO
	1441
[269]	1442	! Potential Temperature
	1443
[1047]	1444	DO ig=1,ngrid
	1445	DO l=1,nlayer
[528]	1446	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
[269]	1447	ENDDO
	1448	ENDDO
	1449
	1450
[226]	1451	c Compute surface stress : (NB: z0 is a common in surfdat.h)
[42]	1452	c DO ig=1,ngrid
[226]	1453	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
[42]	1454	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1455	c ENDDO
	1456
	1457	c Sum of fluxes in solar spectral bands (for output only)
	1458	DO ig=1,ngrid
[226]	1459	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1460	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
[42]	1461	ENDDO
	1462	c ***** TEST ****************************************************
[835]	1463	ztim1 = 999
	1464	DO l=1,nlayer
	1465	DO ig=1,ngrid
	1466	if (pt(ig,l).lt.ztim1) then
	1467	ztim1 = pt(ig,l)
	1468	igmin = ig
	1469	lmin = l
	1470	end if
	1471	ENDDO
	1472	ENDDO
	1473	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
	1474	write(,) 'PHYSIQ: stability WARNING :'
	1475	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1476	& 'ig l =', igmin, lmin
	1477	end if
[42]	1478	c *******************************************************************
	1479
	1480	c ---------------------
	1481	c Outputs to the screen
	1482	c ---------------------
	1483
	1484	IF (lwrite) THEN
	1485	PRINT*,'Global diagnostics for the physics'
	1486	PRINT,'Variables and their increments x and dx/dt dt'
	1487	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1488	WRITE(*,'(2f10.5,2f15.5)')
	1489	s tsurf(igout),zdtsurf(igout)*ptimestep,
[883]	1490	s zplev(igout,1),pdpsrf(igout)*ptimestep
[42]	1491	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1492	WRITE(*,'(i4,6f10.5)') (l,
	1493	s pu(igout,l),pdu(igout,l)*ptimestep,
	1494	s pv(igout,l),pdv(igout,l)*ptimestep,
	1495	s pt(igout,l),pdt(igout,l)*ptimestep,
	1496	s l=1,nlayer)
	1497	ENDIF ! of IF (lwrite)
	1498
[566]	1499	c ----------------------------------------------------------
	1500	c ----------------------------------------------------------
	1501	c INTERPOLATIONS IN THE SURFACE-LAYER
	1502	c ----------------------------------------------------------
	1503	c ----------------------------------------------------------
	1504
[657]	1505	n_out=0 ! number of elements in the z_out array.
[636]	1506	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
	1507	! to 5
	1508	IF (n_out .ne. 0) THEN
[566]	1509
[636]	1510	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
	1511	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
	1512	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
[566]	1513
[636]	1514	z_out(:)=[3.,2.,1.,0.5,0.1]
	1515	u_out(:,:)=0.
	1516	T_out(:,:)=0.
[566]	1517
[636]	1518	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
[566]	1519	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
[657]	1520	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
[566]	1521
	1522	#ifndef MESOSCALE
	1523	IF (ngrid .eq. 1) THEN
	1524	dimout=0
	1525	ELSE
	1526	dimout=2
	1527	ENDIF
	1528	DO n=1,n_out
[636]	1529	write(zstring, '(F8.6)') z_out(n)
	1530	call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
	1531	& 'potential temperature at z_out','K',dimout,T_out(:,n))
[566]	1532	call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
	1533	& 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
[636]	1534	ENDDO
[566]	1535	call WRITEDIAGFI(ngrid,'u_star',
	1536	& 'friction velocity','m/s',dimout,ustar)
	1537	call WRITEDIAGFI(ngrid,'teta_star',
	1538	& 'friction potential temperature','K',dimout,tstar)
	1539	call WRITEDIAGFI(ngrid,'L',
	1540	& 'Monin Obukhov length','m',dimout,L_mo)
	1541	call WRITEDIAGFI(ngrid,'vvv',
	1542	& 'Vertical velocity variance at zout','m',dimout,vvv)
	1543	call WRITEDIAGFI(ngrid,'vhf',
	1544	& 'Vertical heat flux at zout','m',dimout,vhf)
[636]	1545	#else
	1546	T_out1(:)=T_out(:,1)
	1547	u_out1(:)=u_out(:,1)
[566]	1548	#endif
	1549
	1550	ENDIF
	1551
	1552	c ----------------------------------------------------------
	1553	c ----------------------------------------------------------
	1554	c END OF SURFACE LAYER INTERPOLATIONS
	1555	c ----------------------------------------------------------
	1556	c ----------------------------------------------------------
	1557
[42]	1558	IF (ngrid.NE.1) THEN
	1559
[226]	1560	#ifndef MESOSCALE
[42]	1561	c -------------------------------------------------------------------
	1562	c Writing NetCDF file "RESTARTFI" at the end of the run
	1563	c -------------------------------------------------------------------
	1564	c Note: 'restartfi' is stored just before dynamics are stored
	1565	c in 'restart'. Between now and the writting of 'restart',
	1566	c there will have been the itau=itau+1 instruction and
	1567	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1568	c thus we store for time=time+dtvr
	1569
[999]	1570	IF( ((ecritstart.GT.0) .and.
	1571	. (MOD(icount*iphysiq,ecritstart).EQ.0))
	1572	. .or. lastcall ) THEN
	1573
	1574	ztime_fin = pday + ptime + ptimestep/(float(iphysiq)*daysec)
	1575	. - day_ini - time_phys
	1576	print*, pday,ptime,day_ini, time_phys
	1577	write(*,'(A,I7,A,F12.5)')
	1578	. 'PHYSIQ: Ecriture du fichier restartfi ; icount=',
	1579	. icount,' date=',ztime_fin
	1580
	1581
[1047]	1582	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
[999]	1583	. ptimestep,ztime_fin,
	1584	. tsurf,tsoil,co2ice,emis,q2,qsurf)
	1585
[226]	1586	ENDIF
	1587	#endif
[42]	1588
	1589	c -------------------------------------------------------------------
	1590	c Calculation of diagnostic variables written in both stats and
	1591	c diagfi files
	1592	c -------------------------------------------------------------------
	1593
	1594	if (tracer) then
[719]	1595
	1596	if(doubleq) then
	1597	do ig=1,ngrid
	1598	dqdustsurf(ig) =
[756]	1599	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
[719]	1600	dndustsurf(ig) =
[756]	1601	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
	1602	ndust(ig,:) =
	1603	& pq(ig,:,igcm_dust_number)*tauscaling(ig)
	1604	qdust(ig,:) =
	1605	& pq(ig,:,igcm_dust_mass)*tauscaling(ig)
[719]	1606	enddo
	1607	if (scavenging) then
	1608	do ig=1,ngrid
	1609	dqdustsurf(ig) = dqdustsurf(ig) +
[756]	1610	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
[719]	1611	dndustsurf(ig) = dndustsurf(ig) +
[756]	1612	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
	1613	nccn(ig,:) =
	1614	& pq(ig,:,igcm_ccn_number)*tauscaling(ig)
	1615	qccn(ig,:) =
	1616	& pq(ig,:,igcm_ccn_mass)*tauscaling(ig)
[719]	1617	enddo
	1618	endif
	1619	endif
	1620
[42]	1621	if (water) then
[286]	1622	mtot(:)=0
	1623	icetot(:)=0
	1624	rave(:)=0
	1625	tauTES(:)=0
[42]	1626	do ig=1,ngrid
[1047]	1627	do l=1,nlayer
[42]	1628	mtot(ig) = mtot(ig) +
	1629	& zq(ig,l,igcm_h2o_vap) *
[883]	1630	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	1631	icetot(ig) = icetot(ig) +
	1632	& zq(ig,l,igcm_h2o_ice) *
[883]	1633	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	1634	c Computing abs optical depth at 825 cm-1 in each
	1635	c layer to simulate NEW TES retrieval
	1636	Qabsice = min(
	1637	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1638	& )
	1639	opTES(ig,l)= 0.75 * Qabsice *
	1640	& zq(ig,l,igcm_h2o_ice) *
[883]	1641	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	1642	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1643	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1644	enddo
[669]	1645	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
	1646	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
[42]	1647	enddo
[1047]	1648	call watersat(ngrid*nlayer,zt,zplay,zqsat)
[669]	1649	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
[42]	1650
[669]	1651	if (scavenging) then
	1652	Nccntot(:)= 0
	1653	Mccntot(:)= 0
	1654	rave(:)=0
	1655	do ig=1,ngrid
[1047]	1656	do l=1,nlayer
[669]	1657	Nccntot(ig) = Nccntot(ig) +
	1658	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
[883]	1659	& *(zplev(ig,l) - zplev(ig,l+1)) / g
[669]	1660	Mccntot(ig) = Mccntot(ig) +
	1661	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
[883]	1662	& *(zplev(ig,l) - zplev(ig,l+1)) / g
[669]	1663	cccc Column integrated effective ice radius
	1664	cccc is weighted by total ice surface area (BETTER than total ice mass)
	1665	rave(ig) = rave(ig) +
	1666	& tauscaling(ig) *
	1667	& zq(ig,l,igcm_ccn_number) *
[883]	1668	& (zplev(ig,l) - zplev(ig,l+1)) / g *
[669]	1669	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
	1670	enddo
	1671	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
	1672	& /max(pi*rave(ig),1.e-30) ! surface weight
	1673	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1674	enddo
[833]	1675	else ! of if (scavenging)
	1676	rave(:)=0
	1677	do ig=1,ngrid
[1047]	1678	do l=1,nlayer
[833]	1679	rave(ig) = rave(ig) +
	1680	& zq(ig,l,igcm_h2o_ice) *
[883]	1681	& (zplev(ig,l) - zplev(ig,l+1)) / g *
[833]	1682	& rice(ig,l) * (1.+nuice_ref)
	1683	enddo
	1684	rave(ig) = max(rave(ig) /
	1685	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
	1686	enddo
[669]	1687	endif ! of if (scavenging)
	1688
[42]	1689	endif ! of if (water)
	1690	endif ! of if (tracer)
	1691
	1692	c -----------------------------------------------------------------
	1693	c WSTATS: Saving statistics
	1694	c -----------------------------------------------------------------
	1695	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1696	c which can later be used to make the statistic files of the run:
	1697	c "stats") only possible in 3D runs !
	1698
[695]	1699	IF (callstats) THEN
[42]	1700
[695]	1701	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1702	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1703	call wstats(ngrid,"co2ice","CO2 ice cover",
[575]	1704	& "kg.m-2",2,co2ice)
[758]	1705	call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
	1706	& 2,tauref)
[695]	1707	call wstats(ngrid,"fluxsurf_lw",
[575]	1708	& "Thermal IR radiative flux to surface","W.m-2",2,
	1709	& fluxsurf_lw)
[695]	1710	call wstats(ngrid,"fluxsurf_sw",
[575]	1711	& "Solar radiative flux to surface","W.m-2",2,
	1712	& fluxsurf_sw_tot)
[695]	1713	call wstats(ngrid,"fluxtop_lw",
[575]	1714	& "Thermal IR radiative flux to space","W.m-2",2,
	1715	& fluxtop_lw)
[695]	1716	call wstats(ngrid,"fluxtop_sw",
[575]	1717	& "Solar radiative flux to space","W.m-2",2,
	1718	& fluxtop_sw_tot)
[695]	1719	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1720	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1721	call wstats(ngrid,"v","Meridional (North-South) wind",
[575]	1722	& "m.s-1",3,zv)
[705]	1723	call wstats(ngrid,"w","Vertical (down-up) wind",
	1724	& "m.s-1",3,pw)
[695]	1725	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
[883]	1726	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
[758]	1727	call wstats(ngrid,"q2",
	1728	& "Boundary layer eddy kinetic energy",
	1729	& "m2.s-2",3,q2)
	1730	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1731	& emis)
[226]	1732	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1733	c & 2,zstress)
	1734	c call wstats(ngrid,"sw_htrt","sw heat.rate",
	1735	c & "W.m-2",3,zdtsw)
	1736	c call wstats(ngrid,"lw_htrt","lw heat.rate",
	1737	c & "W.m-2",3,zdtlw)
[42]	1738
[758]	1739	if (calltherm) then
	1740	call wstats(ngrid,"zmax_th","Height of thermals",
	1741	& "m",2,zmax_th)
	1742	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
	1743	& "K.m/s",2,hfmax_th)
	1744	call wstats(ngrid,"wstar",
	1745	& "Max vertical velocity in thermals",
	1746	& "m/s",2,wstar)
	1747	endif
	1748
[226]	1749	if (tracer) then
	1750	if (water) then
[1047]	1751	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
	1752	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
[719]	1753	call wstats(ngrid,"vmr_h2ovap",
[520]	1754	& "H2O vapor volume mixing ratio","mol/mol",
	1755	& 3,vmr)
[1047]	1756	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	1757	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
[520]	1758	call wstats(ngrid,"vmr_h2oice",
	1759	& "H2O ice volume mixing ratio","mol/mol",
	1760	& 3,vmr)
[1047]	1761	vmr=zqsat(1:ngrid,1:nlayer)
	1762	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
[669]	1763	call wstats(ngrid,"vmr_h2osat",
	1764	& "saturation volume mixing ratio","mol/mol",
	1765	& 3,vmr)
[226]	1766	call wstats(ngrid,"h2o_ice_s",
	1767	& "surface h2o_ice","kg/m2",
	1768	& 2,qsurf(1,igcm_h2o_ice))
[758]	1769	call wstats(ngrid,'albedo',
	1770	& 'albedo',
[1013]	1771	& '',2,albedo(1,1))
[226]	1772	call wstats(ngrid,"mtot",
	1773	& "total mass of water vapor","kg/m2",
	1774	& 2,mtot)
	1775	call wstats(ngrid,"icetot",
	1776	& "total mass of water ice","kg/m2",
	1777	& 2,icetot)
[520]	1778	call wstats(ngrid,"reffice",
	1779	& "Mean reff","m",
	1780	& 2,rave)
[669]	1781	call wstats(ngrid,"Nccntot",
[520]	1782	& "condensation nuclei","Nbr/m2",
[669]	1783	& 2,Nccntot)
	1784	call wstats(ngrid,"Mccntot",
	1785	& "condensation nuclei mass","kg/m2",
	1786	& 2,Mccntot)
[520]	1787	call wstats(ngrid,"rice",
	1788	& "Ice particle size","m",
	1789	& 3,rice)
[226]	1790	if (.not.activice) then
	1791	call wstats(ngrid,"tauTESap",
	1792	& "tau abs 825 cm-1","",
	1793	& 2,tauTES)
[520]	1794	else
[1047]	1795	call wstats(ngrid,'tauTES',
[520]	1796	& 'tau abs 825 cm-1',
	1797	& '',2,taucloudtes)
[226]	1798	endif
	1799
	1800	endif ! of if (water)
[719]	1801
	1802
	1803	if (dustbin.ne.0) then
	1804
[1047]	1805	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
[719]	1806
	1807	if (doubleq) then
[1047]	1808	c call wstats(ngrid,'qsurf','qsurf',
[719]	1809	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
[1047]	1810	c call wstats(ngrid,'Nsurf','N particles',
[719]	1811	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
[1047]	1812	c call wstats(ngrid,'dqsdev','ddevil lift',
[719]	1813	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
[1047]	1814	c call wstats(ngrid,'dqssed','sedimentation',
[719]	1815	c & 'kg.m-2.s-1',2,zdqssed(1,1))
[1047]	1816	c call wstats(ngrid,'dqsdif','diffusion',
[719]	1817	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
[1047]	1818	call wstats(ngrid,'dqsdust',
[719]	1819	& 'deposited surface dust mass',
	1820	& 'kg.m-2.s-1',2,dqdustsurf)
[1047]	1821	call wstats(ngrid,'dqndust',
[719]	1822	& 'deposited surface dust number',
	1823	& 'number.m-2.s-1',2,dndustsurf)
[1047]	1824	call wstats(ngrid,'reffdust','reffdust',
[719]	1825	& 'm',3,rdust*ref_r0)
[1047]	1826	call wstats(ngrid,'dustq','Dust mass mr',
[756]	1827	& 'kg/kg',3,qdust)
[1047]	1828	call wstats(ngrid,'dustN','Dust number',
[756]	1829	& 'part/kg',3,ndust)
[719]	1830	else
	1831	do iq=1,dustbin
	1832	write(str2(1:2),'(i2.2)') iq
[1047]	1833	call wstats(ngrid,'q'//str2,'mix. ratio',
[719]	1834	& 'kg/kg',3,zq(1,1,iq))
[1047]	1835	call wstats(ngrid,'qsurf'//str2,'qsurf',
[719]	1836	& 'kg.m-2',2,qsurf(1,iq))
	1837	end do
	1838	endif ! (doubleq)
[226]	1839
[719]	1840	if (scavenging) then
[1047]	1841	call wstats(ngrid,'ccnq','CCN mass mr',
[756]	1842	& 'kg/kg',3,qccn)
[1047]	1843	call wstats(ngrid,'ccnN','CCN number',
[756]	1844	& 'part/kg',3,nccn)
[719]	1845	endif ! (scavenging)
	1846
	1847	endif ! (dustbin.ne.0)
	1848
	1849
[226]	1850	if (thermochem.or.photochem) then
	1851	do iq=1,nq
[576]	1852	if (noms(iq) .ne. "dust_mass" .and.
	1853	$ noms(iq) .ne. "dust_number" .and.
	1854	$ noms(iq) .ne. "ccn_mass" .and.
[719]	1855	$ noms(iq) .ne. "ccn_number" .and.
	1856	$ noms(iq) .ne. "h2o_vap" .and.
	1857	$ noms(iq) .ne. "h2o_ice") then
[705]	1858	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	1859	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
	1860	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	1861	& rho(1:ngrid,1:nlayer)n_avog/
	1862	& (1000*mmol(iq))
[576]	1863	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1864	$ "Volume mixing ratio","mol/mol",3,vmr)
[705]	1865	! call wstats(ngrid,"rho_"//trim(noms(iq)),
	1866	! $ "Number density","cm-3",3,rhopart)
	1867	! call writediagfi(ngrid,"rho_"//trim(noms(iq)),
	1868	! $ "Number density","cm-3",3,rhopart)
[576]	1869	if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
	1870	$ (noms(iq).eq."o3")) then
	1871	call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
	1872	$ "Volume mixing ratio","mol/mol",3,vmr)
	1873	end if
	1874	do ig = 1,ngrid
	1875	colden(ig,iq) = 0.
	1876	end do
	1877	do l=1,nlayer
	1878	do ig=1,ngrid
	1879	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
[883]	1880	$ *(zplev(ig,l)-zplev(ig,l+1))
[576]	1881	$ 6.022e22/(mmol(iq)g)
	1882	end do
	1883	end do
	1884	call wstats(ngrid,"c_"//trim(noms(iq)),
	1885	$ "column","mol cm-2",2,colden(1,iq))
	1886	call writediagfi(ngrid,"c_"//trim(noms(iq)),
	1887	$ "column","mol cm-2",2,colden(1,iq))
	1888	end if ! of if (noms(iq) .ne. "dust_mass" ...)
	1889	end do ! of do iq=1,nq
[334]	1890	end if ! of if (thermochem.or.photochem)
[226]	1891
[334]	1892	end if ! of if (tracer)
[226]	1893
	1894	IF(lastcall) THEN
	1895	write (,) "Writing stats..."
	1896	call mkstats(ierr)
	1897	ENDIF
	1898
[42]	1899	ENDIF !if callstats
	1900
	1901	c (Store EOF for Mars Climate database software)
	1902	IF (calleofdump) THEN
	1903	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1904	ENDIF
	1905
[234]	1906
[226]	1907	#ifdef MESOSCALE
[234]	1908	!!!
	1909	!!! OUTPUT FIELDS
	1910	!!!
[81]	1911	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
[668]	1912	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
	1913	TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref)
[667]	1914	IF (tracer) THEN
[82]	1915	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
[608]	1916	icetot(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
[308]	1917	!! JF
	1918	IF (igcm_dust_mass .ne. 0) THEN
[708]	1919	qsurfdust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
[308]	1920	ENDIF
[81]	1921	IF (igcm_h2o_ice .ne. 0) THEN
	1922	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
	1923	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
[1047]	1924	. *mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
[81]	1925	ENDIF
[308]	1926	!! Dust quantity integration along the vertical axe
	1927	dustot(:)=0
[698]	1928	IF (igcm_dust_mass .ne. 0) THEN
[308]	1929	do ig=1,ngrid
[1047]	1930	do l=1,nlayer
[308]	1931	dustot(ig) = dustot(ig) +
	1932	& zq(ig,l,igcm_dust_mass)
[883]	1933	& * (zplev(ig,l) - zplev(ig,l+1)) / g
[308]	1934	enddo
	1935	enddo
[667]	1936	ENDIF
[698]	1937	ENDIF
[550]	1938	!! TAU water ice as seen by TES
	1939	if (activice) tauTES = taucloudtes
[226]	1940	c AUTOMATICALLY GENERATED FROM REGISTRY
[42]	1941	#include "fill_save.inc"
[226]	1942	#else
[528]	1943	#ifndef MESOINI
[42]	1944
[226]	1945	c ==========================================================
	1946	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1947	c any variable for diagnostic (output with period
	1948	c "ecritphy", set in "run.def")
	1949	c ==========================================================
	1950	c WRITEDIAGFI can ALSO be called from any other subroutines
	1951	c for any variables !!
[835]	1952	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1953	c & emis)
	1954	c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
	1955	c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
[226]	1956	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1957	& tsurf)
	1958	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
[769]	1959	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
	1960	& ,"kg.m-2",2,co2ice)
[284]	1961
[719]	1962	call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1963	& "K",2,zt(1,7))
	1964	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1965	& fluxsurf_lw)
	1966	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1967	& fluxsurf_sw_tot)
	1968	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1969	& fluxtop_lw)
	1970	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1971	& fluxtop_sw_tot)
	1972	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1973	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1974	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1975	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1976	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
[226]	1977	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
[520]	1978	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
[758]	1979	call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
[226]	1980	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1981	c & zstress)
[1047]	1982	c call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
[835]	1983	c & 'w.m-2',3,zdtsw)
[1047]	1984	c call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
[835]	1985	c & 'w.m-2',3,zdtlw)
[520]	1986	if (.not.activice) then
[1047]	1987	CALL WRITEDIAGFI(ngrid,'tauTESap',
[520]	1988	& 'tau abs 825 cm-1',
	1989	& '',2,tauTES)
	1990	else
[1047]	1991	CALL WRITEDIAGFI(ngrid,'tauTES',
[520]	1992	& 'tau abs 825 cm-1',
	1993	& '',2,taucloudtes)
	1994	endif
[42]	1995
[528]	1996	#else
	1997	!!! this is to ensure correct initialisation of mesoscale model
	1998	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1999	& tsurf)
	2000	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	2001	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	2002	& co2ice)
	2003	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	2004	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	2005	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
[299]	2006	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	2007	& emis)
	2008	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
	2009	& "K",3,tsoil)
	2010	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
	2011	& "K",3,inertiedat)
	2012	#endif
[42]	2013
[299]	2014
[226]	2015	c ----------------------------------------------------------
	2016	c Outputs of the CO2 cycle
	2017	c ----------------------------------------------------------
[42]	2018
[226]	2019	if (tracer.and.(igcm_co2.ne.0)) then
	2020	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	2021	! & "kg/kg",2,zq(1,1,igcm_co2))
[719]	2022	call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	2023	& "kg/kg",3,zq(1,1,igcm_co2))
[226]	2024
	2025	! Compute co2 column
[286]	2026	co2col(:)=0
[1047]	2027	do l=1,nlayer
[226]	2028	do ig=1,ngrid
	2029	co2col(ig)=co2col(ig)+
[883]	2030	& zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
[226]	2031	enddo
	2032	enddo
	2033	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	2034	& co2col)
	2035	endif ! of if (tracer.and.(igcm_co2.ne.0))
[42]	2036
[226]	2037	c ----------------------------------------------------------
	2038	c Outputs of the water cycle
	2039	c ----------------------------------------------------------
	2040	if (tracer) then
	2041	if (water) then
	2042
[299]	2043	#ifdef MESOINI
	2044	!!!! waterice = q01, voir readmeteo.F90
[1047]	2045	call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
[299]	2046	& 'kg/kg',3,
[1047]	2047	& zq(1:ngrid,1:nlayer,igcm_h2o_ice))
[299]	2048	!!!! watervapor = q02, voir readmeteo.F90
[1047]	2049	call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
[299]	2050	& 'kg/kg',3,
[1047]	2051	& zq(1:ngrid,1:nlayer,igcm_h2o_vap))
[299]	2052	!!!! surface waterice qsurf02 (voir readmeteo)
[1047]	2053	call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
[299]	2054	& 'kg.m-2',2,
[1047]	2055	& qsurf(1:ngrid,igcm_h2o_ice))
[299]	2056	#endif
	2057
[1047]	2058	CALL WRITEDIAGFI(ngrid,'mtot',
[226]	2059	& 'total mass of water vapor',
	2060	& 'kg/m2',2,mtot)
[1047]	2061	CALL WRITEDIAGFI(ngrid,'icetot',
[226]	2062	& 'total mass of water ice',
	2063	& 'kg/m2',2,icetot)
[1047]	2064	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	2065	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
	2066	call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
[719]	2067	& 'mol/mol',3,vmr)
[1047]	2068	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
	2069	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
	2070	call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
[719]	2071	& 'mol/mol',3,vmr)
[1047]	2072	CALL WRITEDIAGFI(ngrid,'reffice',
[669]	2073	& 'Mean reff',
	2074	& 'm',2,rave)
[719]	2075	CALL WRITEDIAGFI(ngrid,"Nccntot",
[669]	2076	& "condensation nuclei","Nbr/m2",
	2077	& 2,Nccntot)
[719]	2078	CALL WRITEDIAGFI(ngrid,"Mccntot",
	2079	& "mass condensation nuclei","kg/m2",
	2080	& 2,Mccntot)
[1047]	2081	call WRITEDIAGFI(ngrid,'rice','Ice particle size',
[719]	2082	& 'm',3,rice)
[1047]	2083	call WRITEDIAGFI(ngrid,'h2o_ice_s',
[226]	2084	& 'surface h2o_ice',
	2085	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
[1047]	2086	CALL WRITEDIAGFI(ngrid,'albedo',
[833]	2087	& 'albedo',
[1013]	2088	& '',2,albedo(1,1))
[833]	2089	if (tifeedback) then
[1047]	2090	call WRITEDIAGSOIL(ngrid,"soiltemp",
[833]	2091	& "Soil temperature","K",
	2092	& 3,tsoil)
[1047]	2093	call WRITEDIAGSOIL(ngrid,'soilti',
[833]	2094	& 'Soil Thermal Inertia',
	2095	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
	2096	endif
[226]	2097	endif !(water)
	2098
	2099
	2100	if (water.and..not.photochem) then
	2101	iq=nq
	2102	c write(str2(1:2),'(i2.2)') iq
[1047]	2103	c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
[226]	2104	c & 'kg.m-2',2,zdqscloud(1,iq))
[1047]	2105	c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
[226]	2106	c & 'kg/kg',3,zdqchim(1,1,iq))
[1047]	2107	c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
[226]	2108	c & 'kg/kg',3,zdqdif(1,1,iq))
[1047]	2109	c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
[226]	2110	c & 'kg/kg',3,zdqadj(1,1,iq))
[1047]	2111	c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
[226]	2112	c & 'kg/kg',3,zdqc(1,1,iq))
	2113	endif !(water.and..not.photochem)
	2114	endif
	2115
	2116	c ----------------------------------------------------------
	2117	c Outputs of the dust cycle
	2118	c ----------------------------------------------------------
	2119
[1047]	2120	call WRITEDIAGFI(ngrid,'tauref',
[758]	2121	& 'Dust ref opt depth','NU',2,tauref)
[226]	2122
	2123	if (tracer.and.(dustbin.ne.0)) then
[1047]	2124	call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
[226]	2125	if (doubleq) then
[1047]	2126	c call WRITEDIAGFI(ngrid,'qsurf','qsurf',
[411]	2127	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
[1047]	2128	c call WRITEDIAGFI(ngrid,'Nsurf','N particles',
[411]	2129	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
[1047]	2130	c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
[226]	2131	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
[1047]	2132	c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
[411]	2133	c & 'kg.m-2.s-1',2,zdqssed(1,1))
[1047]	2134	c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
[411]	2135	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
[1047]	2136	call WRITEDIAGFI(ngrid,'dqsdust',
[719]	2137	& 'deposited surface dust mass',
	2138	& 'kg.m-2.s-1',2,dqdustsurf)
[1047]	2139	call WRITEDIAGFI(ngrid,'dqndust',
[719]	2140	& 'deposited surface dust number',
	2141	& 'number.m-2.s-1',2,dndustsurf)
[1047]	2142	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
[719]	2143	& 'm',3,rdust*ref_r0)
[1047]	2144	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
[756]	2145	& 'kg/kg',3,qdust)
[1047]	2146	call WRITEDIAGFI(ngrid,'dustN','Dust number',
[756]	2147	& 'part/kg',3,ndust)
[299]	2148	#ifdef MESOINI
[1032]	2149	! !!! to initialize mesoscale we need scaled variables
	2150	! !!! because this must correspond to starting point for tracers
[1047]	2151	! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
	2152	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
	2153	! call WRITEDIAGFI(ngrid,'dustN','Dust number',
	2154	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
	2155	! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
	2156	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
	2157	! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
	2158	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
	2159	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
[1032]	2160	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
[1047]	2161	call WRITEDIAGFI(ngrid,'dustN','Dust number',
[1032]	2162	& 'part/kg',3,pq(1,1,igcm_dust_number))
[1047]	2163	call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
[1032]	2164	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
[1047]	2165	call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
[1032]	2166	& 'part/kg',3,pq(1,1,igcm_ccn_number))
[299]	2167	#endif
[226]	2168	else
	2169	do iq=1,dustbin
	2170	write(str2(1:2),'(i2.2)') iq
[1047]	2171	call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
[226]	2172	& 'kg/kg',3,zq(1,1,iq))
[1047]	2173	call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
[226]	2174	& 'kg.m-2',2,qsurf(1,iq))
	2175	end do
	2176	endif ! (doubleq)
[358]	2177
	2178	if (scavenging) then
[1047]	2179	call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
[756]	2180	& 'kg/kg',3,qccn)
[1047]	2181	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
[756]	2182	& 'part/kg',3,nccn)
[358]	2183	endif ! (scavenging)
	2184
[226]	2185	c if (submicron) then
[1047]	2186	c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
[226]	2187	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
	2188	c endif ! (submicron)
	2189	end if ! (tracer.and.(dustbin.ne.0))
	2190
[705]	2191
[226]	2192	c ----------------------------------------------------------
[705]	2193	c Thermospheric outputs
[267]	2194	c ----------------------------------------------------------
[705]	2195
	2196	if(callthermos) then
	2197
[1047]	2198	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
[705]	2199	$ 3,zdtnlte)
[1047]	2200	call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
[705]	2201	$ 3,zdteuv)
[1047]	2202	call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
[705]	2203	$ 3,zdtconduc)
[1047]	2204	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
[705]	2205	$ 3,zdtnirco2)
	2206
	2207	endif !(callthermos)
	2208
	2209	c ----------------------------------------------------------
	2210	c ----------------------------------------------------------
[267]	2211	c PBL OUTPUS
	2212	c ----------------------------------------------------------
	2213	c ----------------------------------------------------------
	2214
	2215	c ----------------------------------------------------------
[226]	2216	c Outputs of thermals
	2217	c ----------------------------------------------------------
	2218	if (calltherm) then
	2219
	2220	! call WRITEDIAGFI(ngrid,'dtke',
	2221	! & 'tendance tke thermiques','m2/s2',
	2222	! & 3,dtke_th)
	2223	! call WRITEDIAGFI(ngrid,'d_u_ajs',
	2224	! & 'tendance u thermiques','m/s',
	2225	! & 3,pdu_th*ptimestep)
	2226	! call WRITEDIAGFI(ngrid,'d_v_ajs',
	2227	! & 'tendance v thermiques','m/s',
	2228	! & 3,pdv_th*ptimestep)
	2229	! if (tracer) then
	2230	! if (nq .eq. 2) then
	2231	! call WRITEDIAGFI(ngrid,'deltaq_th',
	2232	! & 'delta q thermiques','kg/kg',
	2233	! & 3,ptimestep*pdq_th(:,:,2))
	2234	! endif
	2235	! endif
	2236
[1047]	2237	call WRITEDIAGFI(ngrid,'zmax_th',
[523]	2238	& 'hauteur du thermique','m',
	2239	& 2,zmax_th)
[1047]	2240	call WRITEDIAGFI(ngrid,'hfmax_th',
[226]	2241	& 'maximum TH heat flux','K.m/s',
	2242	& 2,hfmax_th)
[1047]	2243	call WRITEDIAGFI(ngrid,'wstar',
[226]	2244	& 'maximum TH vertical velocity','m/s',
[528]	2245	& 2,wstar)
[226]	2246
	2247	endif
	2248
	2249	c ----------------------------------------------------------
[267]	2250	c ----------------------------------------------------------
	2251	c END OF PBL OUTPUS
	2252	c ----------------------------------------------------------
	2253	c ----------------------------------------------------------
	2254
	2255
	2256	c ----------------------------------------------------------
[226]	2257	c Output in netcdf file "diagsoil.nc" for subterranean
	2258	c variables (output every "ecritphy", as for writediagfi)
	2259	c ----------------------------------------------------------
	2260
	2261	! Write soil temperature
	2262	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	2263	! & 3,tsoil)
	2264	! Write surface temperature
	2265	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	2266	! & 2,tsurf)
	2267
	2268	c ==========================================================
	2269	c END OF WRITEDIAGFI
	2270	c ==========================================================
	2271	#endif
	2272
[42]	2273	ELSE ! if(ngrid.eq.1)
	2274
[627]	2275	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
	2276	& zls*180./pi,odpref,tauref
[42]	2277	c ----------------------------------------------------------------------
	2278	c Output in grads file "g1d" (ONLY when using testphys1d)
	2279	c (output at every X physical timestep)
	2280	c ----------------------------------------------------------------------
	2281	c
	2282	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	2283	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	2284	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	2285
	2286	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	2287	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	2288	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	2289	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	2290
[226]	2291	! THERMALS STUFF 1D
	2292	if(calltherm) then
	2293
[1047]	2294	call WRITEDIAGFI(ngrid,'lmax_th',
[226]	2295	& 'hauteur du thermique','point',
	2296	& 0,lmax_th_out)
[1047]	2297	call WRITEDIAGFI(ngrid,'zmax_th',
[528]	2298	& 'hauteur du thermique','m',
	2299	& 0,zmax_th)
[1047]	2300	call WRITEDIAGFI(ngrid,'hfmax_th',
[226]	2301	& 'maximum TH heat flux','K.m/s',
	2302	& 0,hfmax_th)
[1047]	2303	call WRITEDIAGFI(ngrid,'wstar',
[226]	2304	& 'maximum TH vertical velocity','m/s',
[528]	2305	& 0,wstar)
[226]	2306
	2307	co2col(:)=0.
	2308	if (tracer) then
[1047]	2309	do l=1,nlayer
[226]	2310	do ig=1,ngrid
	2311	co2col(ig)=co2col(ig)+
[883]	2312	& zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
[226]	2313	enddo
	2314	enddo
	2315
	2316	end if
	2317	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
	2318	& ,'kg/m-2',0,co2col)
[722]	2319	endif ! of if (calltherm)
	2320
[226]	2321	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
	2322	& ,'m/s',1,pw)
	2323	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
	2324	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
	2325	& tsurf)
[277]	2326	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
	2327	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
[226]	2328
	2329	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
	2330	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
[358]	2331	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
[544]	2332	! call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
	2333	! & "K.s-1",1,dtrad/zpopsk)
[1047]	2334	! call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
[544]	2335	! & 'w.m-2',1,zdtsw/zpopsk)
[1047]	2336	! call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
[544]	2337	! & 'w.m-2',1,zdtlw/zpopsk)
[769]	2338	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
	2339	& ,"kg.m-2",0,co2ice)
[544]	2340
[226]	2341	! or output in diagfi.nc (for testphys1d)
[1047]	2342	call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
	2343	call WRITEDIAGFI(ngrid,'temp','Temperature',
[226]	2344	& 'K',1,zt)
[358]	2345
[226]	2346	if(tracer) then
	2347	c CALL writeg1d(ngrid,1,tau,'tau','SI')
	2348	do iq=1,nq
	2349	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
[1047]	2350	call WRITEDIAGFI(ngrid,trim(noms(iq)),
[226]	2351	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	2352	end do
[358]	2353	if (doubleq) then
[1047]	2354	call WRITEDIAGFI(ngrid,'rdust','rdust',
[358]	2355	& 'm',1,rdust)
	2356	endif
[833]	2357	if (water.AND.tifeedback) then
[1047]	2358	call WRITEDIAGFI(ngrid,"soiltemp",
[833]	2359	& "Soil temperature","K",
	2360	& 1,tsoil)
[1047]	2361	call WRITEDIAGFI(ngrid,'soilti',
[833]	2362	& 'Soil Thermal Inertia',
	2363	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
	2364	endif
[226]	2365	end if
[358]	2366
[520]	2367	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
[358]	2368	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[520]	2369	IF (water) THEN
[722]	2370
	2371	if (.not.activice) then
	2372
	2373	tauTES=0
[1047]	2374	do l=1,nlayer
[722]	2375	Qabsice = min(
	2376	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
	2377	& )
	2378	opTES(1,l)= 0.75 * Qabsice *
	2379	& zq(1,l,igcm_h2o_ice) *
[883]	2380	& (zplev(1,l) - zplev(1,l+1)) / g
[722]	2381	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
	2382	tauTES=tauTES+ opTES(1,l)
	2383	enddo
[1047]	2384	CALL WRITEDIAGFI(ngrid,'tauTESap',
[358]	2385	& 'tau abs 825 cm-1',
[520]	2386	& '',0,tauTES)
[722]	2387	else
	2388
[1047]	2389	CALL WRITEDIAGFI(ngrid,'tauTES',
[520]	2390	& 'tau abs 825 cm-1',
	2391	& '',0,taucloudtes)
[722]	2392	endif
[520]	2393
	2394	mtot = 0
	2395	icetot = 0
	2396	h2otot = qsurf(1,igcm_h2o_ice)
[722]	2397
[358]	2398	do l=1,nlayer
[520]	2399	mtot = mtot + zq(1,l,igcm_h2o_vap)
[883]	2400	& * (zplev(1,l) - zplev(1,l+1)) / g
[520]	2401	icetot = icetot + zq(1,l,igcm_h2o_ice)
[883]	2402	& * (zplev(1,l) - zplev(1,l+1)) / g
[722]	2403	end do
	2404	h2otot = h2otot+mtot+icetot
	2405
[1047]	2406	CALL WRITEDIAGFI(ngrid,'h2otot',
[722]	2407	& 'h2otot',
	2408	& 'kg/m2',0,h2otot)
[1047]	2409	CALL WRITEDIAGFI(ngrid,'mtot',
[722]	2410	& 'mtot',
	2411	& 'kg/m2',0,mtot)
[1047]	2412	CALL WRITEDIAGFI(ngrid,'icetot',
[722]	2413	& 'icetot',
	2414	& 'kg/m2',0,icetot)
	2415
	2416	if (scavenging) then
	2417
	2418	rave = 0
	2419	do l=1,nlayer
[669]	2420	cccc Column integrated effective ice radius
	2421	cccc is weighted by total ice surface area (BETTER)
[722]	2422	rave = rave + tauscaling(1) *
	2423	& zq(1,l,igcm_ccn_number) *
[883]	2424	& (zplev(1,l) - zplev(1,l+1)) / g *
[722]	2425	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
	2426	enddo
	2427	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
	2428
	2429	Nccntot= 0
[1047]	2430	call watersat(ngrid*nlayer,zt,zplay,zqsat)
	2431	do l=1,nlayer
[669]	2432	Nccntot = Nccntot +
[520]	2433	& zq(1,l,igcm_ccn_number)*tauscaling(1)
[883]	2434	& *(zplev(1,l) - zplev(1,l+1)) / g
[520]	2435	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
	2436	satu(1,l) = (max(satu(1,l),float(1))-1)
	2437	! & * zq(1,l,igcm_h2o_vap) *
[883]	2438	! & (zplev(1,l) - zplev(1,l+1)) / g
[520]	2439	enddo
[722]	2440	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
	2441	& satu)
[1047]	2442	CALL WRITEDIAGFI(ngrid,'Nccntot',
[669]	2443	& 'Nccntot',
	2444	& 'nbr/m2',0,Nccntot)
[411]	2445
[1047]	2446	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
[769]	2447	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
[1047]	2448	call WRITEDIAGFI(ngrid,'zdqsed_dustN'
[769]	2449	&,'sedimentation N','Nbr.m-2.s-1',1,
	2450	& zdqsed(1,:,igcm_dust_number))
[722]	2451
	2452	else ! of if (scavenging)
	2453
	2454	cccc Column integrated effective ice radius
	2455	cccc is weighted by total ice mass (LESS GOOD)
	2456	rave = 0
	2457	do l=1,nlayer
	2458	rave = rave + zq(1,l,igcm_h2o_ice)
[883]	2459	& * (zplev(1,l) - zplev(1,l+1)) / g
[722]	2460	& * rice(1,l) * (1.+nuice_ref)
	2461	enddo
	2462	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
	2463	endif ! of if (scavenging)
	2464
[1047]	2465	CALL WRITEDIAGFI(ngrid,'reffice',
[722]	2466	& 'reffice',
	2467	& 'm',0,rave)
	2468
	2469	do iq=1,nq
[1047]	2470	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
[722]	2471	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
	2472	end do
[358]	2473
[1047]	2474	call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
[633]	2475	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
[1047]	2476	call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
[633]	2477	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
[1047]	2478	call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
[633]	2479	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
	2480	& +zdqcloud(1,:,igcm_h2o_vap))
[411]	2481
[520]	2482	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
[411]	2483	& rice)
[520]	2484	ENDIF ! of IF (water)
[411]	2485
[358]	2486	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2487	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2488
	2489
[883]	2490	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
[226]	2491
	2492	do l=2,nlayer-1
	2493	tmean=zt(1,l)
	2494	if(zt(1,l).ne.zt(1,l-1))
	2495	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	2496	zlocal(l)= zlocal(l-1)
[883]	2497	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
[226]	2498	enddo
	2499	zlocal(nlayer)= zlocal(nlayer-1)-
[883]	2500	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
[226]	2501	& rnew(1,nlayer)*tmean/g
	2502
[42]	2503	END IF ! if(ngrid.ne.1)
	2504
	2505	icount=icount+1
	2506	RETURN
	2507	END

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

Download in other formats:

Original Format