Context Navigation

physiq.F @ 708

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

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

Download in other formats:

Original Format