Context Navigation

physiq.F @ 1047

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

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

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

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

Download in other formats:

Original Format