Context Navigation

physiq.F @ 756

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

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

Download in other formats:

Original Format