Context Navigation

physiq.F @ 705

Last change on this file since 705 was 705, checked in by emillour, 13 years ago

Mars GCM:

Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

FGG

File size: 81.9 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
[633]	1006	pdq(1:ngrid,1:nlayer,1:nq) =
	1007	& pdq(1:ngrid,1:nlayer,1:nq) +
	1008	& zdqcloud(1:ngrid,1:nlayer,1:nq)
	1009
	1010	! We need to check that we have Nccn & Ndust > 0
	1011	! This is due to single precision rounding problems
	1012	if (scavenging) then
[411]	1013
[633]	1014	where (pq(:,:,igcm_ccn_mass) +
	1015	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
	1016	pdq(:,:,igcm_ccn_mass) =
	1017	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
	1018	pdq(:,:,igcm_ccn_number) =
	1019	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1020	end where
	1021	where (pq(:,:,igcm_ccn_number) +
	1022	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
	1023	pdq(:,:,igcm_ccn_mass) =
	1024	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
	1025	pdq(:,:,igcm_ccn_number) =
	1026	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1027	end where
	1028	where (pq(:,:,igcm_dust_mass) +
	1029	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
	1030	pdq(:,:,igcm_dust_mass) =
	1031	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
	1032	pdq(:,:,igcm_dust_number) =
	1033	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1034	end where
	1035	where (pq(:,:,igcm_dust_number) +
	1036	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
	1037	pdq(:,:,igcm_dust_mass) =
	1038	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
	1039	pdq(:,:,igcm_dust_number) =
	1040	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1041	end where
	1042
	1043	endif ! of if scavenging
	1044
	1045
[42]	1046	END IF ! of IF (water)
	1047
[556]	1048	c 7b. Aerosol particles
[42]	1049	c -------------------
	1050
	1051	c ----------
	1052	c Dust devil :
	1053	c ----------
	1054	IF(callddevil) then
	1055	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
	1056	& zdqdev,zdqsdev)
	1057
	1058	if (dustbin.ge.1) then
	1059	do iq=1,nq
	1060	DO l=1,nlayer
	1061	DO ig=1,ngrid
	1062	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	1063	ENDDO
	1064	ENDDO
	1065	enddo
	1066	do iq=1,nq
	1067	DO ig=1,ngrid
	1068	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
	1069	ENDDO
	1070	enddo
	1071	endif ! of if (dustbin.ge.1)
	1072
	1073	END IF ! of IF (callddevil)
	1074
	1075	c -------------
	1076	c Sedimentation : acts also on water ice
	1077	c -------------
	1078	IF (sedimentation) THEN
	1079	!call zerophys(ngridnlayernq, zdqsed)
	1080	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	1081	!call zerophys(ngrid*nq, zdqssed)
	1082	zdqssed(1:ngrid,1:nq)=0
	1083
	1084	call callsedim(ngrid,nlayer, ptimestep,
[411]	1085	& pplev,zzlev, zzlay, pt, rdust, rice,
[358]	1086	& rsedcloud,rhocloud,
[411]	1087	& pq, pdq, zdqsed, zdqssed,nq,
	1088	& tau,tauscaling)
[520]	1089
	1090
[42]	1091	DO iq=1, nq
	1092	DO l=1,nlayer
	1093	DO ig=1,ngrid
	1094	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1095	ENDDO
	1096	ENDDO
	1097	ENDDO
	1098	DO iq=1, nq
	1099	DO ig=1,ngrid
	1100	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1101	ENDDO
	1102	ENDDO
	1103	END IF ! of IF (sedimentation)
[520]	1104
[556]	1105	c
	1106	c 7c. Chemical species
	1107	c ------------------
[42]	1108
[556]	1109	#ifndef MESOSCALE
	1110	c --------------
	1111	c photochemistry :
	1112	c --------------
	1113	IF (photochem .or. thermochem) then
[411]	1114
[556]	1115	! dust and ice surface area
	1116	call surfacearea(ngrid, nlayer, ptimestep, pplay, zzlay,
	1117	$ pt, pq, pdq, nq,
	1118	$ rdust, rice, tau, tauscaling,
	1119	$ surfdust, surfice)
	1120	! call photochemistry
	1121	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
	1122	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
	1123	$ zdqcloud,zdqscloud,tauref,co2ice,
	1124	$ pu,pdu,pv,pdv,surfdust,surfice)
	1125
	1126	! increment values of tracers:
	1127	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1128	! tracers is zero anyways
	1129	DO l=1,nlayer
	1130	DO ig=1,ngrid
	1131	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	1132	ENDDO
	1133	ENDDO
	1134	ENDDO ! of DO iq=1,nq
	1135
	1136	! add condensation tendency for H2O2
	1137	if (igcm_h2o2.ne.0) then
	1138	DO l=1,nlayer
	1139	DO ig=1,ngrid
	1140	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	1141	& +zdqcloud(ig,l,igcm_h2o2)
	1142	ENDDO
	1143	ENDDO
	1144	endif
	1145
	1146	! increment surface values of tracers:
	1147	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1148	! tracers is zero anyways
	1149	DO ig=1,ngrid
	1150	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	1151	ENDDO
	1152	ENDDO ! of DO iq=1,nq
	1153
	1154	! add condensation tendency for H2O2
	1155	if (igcm_h2o2.ne.0) then
	1156	DO ig=1,ngrid
	1157	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	1158	& +zdqscloud(ig,igcm_h2o2)
	1159	ENDDO
	1160	endif
	1161
	1162	END IF ! of IF (photochem.or.thermochem)
	1163	#endif
	1164
[42]	1165	c 7d. Updates
	1166	c ---------
	1167
	1168	DO iq=1, nq
	1169	DO ig=1,ngrid
	1170
	1171	c ---------------------------------
	1172	c Updating tracer budget on surface
	1173	c ---------------------------------
	1174	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1175
	1176	ENDDO ! (ig)
	1177	ENDDO ! (iq)
	1178
	1179	endif ! of if (tracer)
	1180
[226]	1181	#ifndef MESOSCALE
	1182	c-----------------------------------------------------------------------
	1183	c 8. THERMOSPHERE CALCULATION
	1184	c-----------------------------------------------------------------------
[42]	1185
[226]	1186	if (callthermos) then
	1187	call thermosphere(pplev,pplay,dist_sol,
	1188	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	1189	& pt,pq,pu,pv,pdt,pdq,
	1190	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	1191
	1192	DO l=1,nlayer
	1193	DO ig=1,ngrid
	1194	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	1195	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	1196	& +zdteuv(ig,l)
	1197	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	1198	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	1199	DO iq=1, nq
	1200	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	1201	ENDDO
	1202	ENDDO
	1203	ENDDO
	1204
	1205	endif ! of if (callthermos)
	1206	#endif
	1207
[42]	1208	c-----------------------------------------------------------------------
	1209	c 9. Surface and sub-surface soil temperature
	1210	c-----------------------------------------------------------------------
	1211	c
	1212	c
	1213	c 9.1 Increment Surface temperature:
	1214	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1215
	1216	DO ig=1,ngrid
	1217	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1218	ENDDO
	1219
	1220	c Prescribe a cold trap at south pole (except at high obliquity !!)
	1221	c Temperature at the surface is set there to be the temperature
	1222	c corresponding to equilibrium temperature between phases of CO2
	1223
[411]	1224
[42]	1225	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
[226]	1226	#ifndef MESOSCALE
	1227	if (caps.and.(obliquit.lt.27.)) then
	1228	! NB: Updated surface pressure, at grid point 'ngrid', is
	1229	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
	1230	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	1231	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	1232	endif
	1233	#endif
[42]	1234	c -------------------------------------------------------------
[283]	1235	c Change of surface albedo in case of ground frost
[42]	1236	c everywhere except on the north permanent cap and in regions
	1237	c covered by dry ice.
	1238	c ALWAYS PLACE these lines after newcondens !!!
	1239	c -------------------------------------------------------------
	1240	do ig=1,ngrid
	1241	if ((co2ice(ig).eq.0).and.
[283]	1242	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
	1243	albedo(ig,1) = albedo_h2o_ice
	1244	albedo(ig,2) = albedo_h2o_ice
	1245	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
	1246	c write(,) "physiq.F frost :"
	1247	c & ,lati(ig)180./pi, long(ig)180./pi
[42]	1248	endif
	1249	enddo ! of do ig=1,ngrid
	1250	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
	1251
[528]	1252	c
[42]	1253	c 9.2 Compute soil temperatures and subsurface heat flux:
	1254	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1255	IF (callsoil) THEN
	1256	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
	1257	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1258	ENDIF
[411]	1259
[42]	1260
	1261	c-----------------------------------------------------------------------
	1262	c 10. Write output files
	1263	c ----------------------
	1264
[674]	1265	c Save variables for eventual restart in MMM and LES
	1266	#ifdef MESOSCALE
	1267	#include "meso_inc/meso_inc_save_restart.F"
	1268	#endif
	1269
[42]	1270	c -------------------------------
	1271	c Dynamical fields incrementation
	1272	c -------------------------------
	1273	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1274	! temperature, zonal and meridional wind
	1275	DO l=1,nlayer
	1276	DO ig=1,ngrid
	1277	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1278	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1279	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1280	ENDDO
	1281	ENDDO
	1282
	1283	! tracers
	1284	DO iq=1, nq
	1285	DO l=1,nlayer
	1286	DO ig=1,ngrid
	1287	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1288	ENDDO
	1289	ENDDO
	1290	ENDDO
	1291
	1292	! surface pressure
	1293	DO ig=1,ngrid
	1294	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1295	ENDDO
	1296
	1297	! pressure
	1298	DO l=1,nlayer
	1299	DO ig=1,ngrid
	1300	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1301	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1302	ENDDO
	1303	ENDDO
	1304
	1305	! Density
	1306	DO l=1,nlayer
	1307	DO ig=1,ngrid
	1308	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1309	ENDDO
	1310	ENDDO
	1311
[269]	1312	! Potential Temperature
	1313
	1314	DO ig=1,ngridmx
	1315	DO l=1,nlayermx
[528]	1316	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
[269]	1317	ENDDO
	1318	ENDDO
	1319
	1320
[226]	1321	c Compute surface stress : (NB: z0 is a common in surfdat.h)
[42]	1322	c DO ig=1,ngrid
[226]	1323	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
[42]	1324	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1325	c ENDDO
	1326
	1327	c Sum of fluxes in solar spectral bands (for output only)
	1328	DO ig=1,ngrid
[226]	1329	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1330	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
[42]	1331	ENDDO
	1332	c ***** TEST ****************************************************
	1333	ztim1 = 999
	1334	DO l=1,nlayer
	1335	DO ig=1,ngrid
	1336	if (pt(ig,l).lt.ztim1) then
	1337	ztim1 = pt(ig,l)
	1338	igmin = ig
	1339	lmin = l
	1340	end if
	1341	ENDDO
	1342	ENDDO
[234]	1343	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
[42]	1344	write(,) 'PHYSIQ: stability WARNING :'
	1345	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1346	& 'ig l =', igmin, lmin
	1347	end if
	1348	c *******************************************************************
	1349
	1350	c ---------------------
	1351	c Outputs to the screen
	1352	c ---------------------
	1353
	1354	IF (lwrite) THEN
	1355	PRINT*,'Global diagnostics for the physics'
	1356	PRINT,'Variables and their increments x and dx/dt dt'
	1357	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1358	WRITE(*,'(2f10.5,2f15.5)')
	1359	s tsurf(igout),zdtsurf(igout)*ptimestep,
	1360	s pplev(igout,1),pdpsrf(igout)*ptimestep
	1361	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1362	WRITE(*,'(i4,6f10.5)') (l,
	1363	s pu(igout,l),pdu(igout,l)*ptimestep,
	1364	s pv(igout,l),pdv(igout,l)*ptimestep,
	1365	s pt(igout,l),pdt(igout,l)*ptimestep,
	1366	s l=1,nlayer)
	1367	ENDIF ! of IF (lwrite)
	1368
[566]	1369	c ----------------------------------------------------------
	1370	c ----------------------------------------------------------
	1371	c INTERPOLATIONS IN THE SURFACE-LAYER
	1372	c ----------------------------------------------------------
	1373	c ----------------------------------------------------------
	1374
[657]	1375	n_out=0 ! number of elements in the z_out array.
[636]	1376	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
	1377	! to 5
	1378	IF (n_out .ne. 0) THEN
[566]	1379
[636]	1380	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
	1381	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
	1382	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
[566]	1383
[636]	1384	z_out(:)=[3.,2.,1.,0.5,0.1]
	1385	u_out(:,:)=0.
	1386	T_out(:,:)=0.
[566]	1387
[636]	1388	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
[566]	1389	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
[657]	1390	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
[566]	1391
	1392	#ifndef MESOSCALE
	1393	IF (ngrid .eq. 1) THEN
	1394	dimout=0
	1395	ELSE
	1396	dimout=2
	1397	ENDIF
	1398	DO n=1,n_out
[636]	1399	write(zstring, '(F8.6)') z_out(n)
	1400	call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
	1401	& 'potential temperature at z_out','K',dimout,T_out(:,n))
[566]	1402	call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
	1403	& 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
[636]	1404	ENDDO
[566]	1405	call WRITEDIAGFI(ngrid,'u_star',
	1406	& 'friction velocity','m/s',dimout,ustar)
	1407	call WRITEDIAGFI(ngrid,'teta_star',
	1408	& 'friction potential temperature','K',dimout,tstar)
	1409	call WRITEDIAGFI(ngrid,'L',
	1410	& 'Monin Obukhov length','m',dimout,L_mo)
	1411	call WRITEDIAGFI(ngrid,'vvv',
	1412	& 'Vertical velocity variance at zout','m',dimout,vvv)
	1413	call WRITEDIAGFI(ngrid,'vhf',
	1414	& 'Vertical heat flux at zout','m',dimout,vhf)
[636]	1415	#else
	1416	T_out1(:)=T_out(:,1)
	1417	u_out1(:)=u_out(:,1)
[566]	1418	#endif
	1419
	1420	ENDIF
	1421
	1422	c ----------------------------------------------------------
	1423	c ----------------------------------------------------------
	1424	c END OF SURFACE LAYER INTERPOLATIONS
	1425	c ----------------------------------------------------------
	1426	c ----------------------------------------------------------
	1427
[42]	1428	IF (ngrid.NE.1) THEN
	1429
[226]	1430	#ifndef MESOSCALE
[42]	1431	c -------------------------------------------------------------------
	1432	c Writing NetCDF file "RESTARTFI" at the end of the run
	1433	c -------------------------------------------------------------------
	1434	c Note: 'restartfi' is stored just before dynamics are stored
	1435	c in 'restart'. Between now and the writting of 'restart',
	1436	c there will have been the itau=itau+1 instruction and
	1437	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1438	c thus we store for time=time+dtvr
	1439
[226]	1440	IF(lastcall) THEN
	1441	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1442	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1443	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1444	. ptimestep,pday,
	1445	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
	1446	. area,albedodat,inertiedat,zmea,zstd,zsig,
	1447	. zgam,zthe)
	1448	ENDIF
	1449	#endif
[42]	1450
	1451	c -------------------------------------------------------------------
	1452	c Calculation of diagnostic variables written in both stats and
	1453	c diagfi files
	1454	c -------------------------------------------------------------------
	1455
	1456	if (tracer) then
	1457	if (water) then
	1458
[286]	1459	mtot(:)=0
	1460	icetot(:)=0
	1461	rave(:)=0
	1462	tauTES(:)=0
[42]	1463	do ig=1,ngrid
	1464	do l=1,nlayermx
	1465	mtot(ig) = mtot(ig) +
	1466	& zq(ig,l,igcm_h2o_vap) *
	1467	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1468	icetot(ig) = icetot(ig) +
	1469	& zq(ig,l,igcm_h2o_ice) *
	1470	& (pplev(ig,l) - pplev(ig,l+1)) / g
[669]	1471	cccc Column integrated effective ice radius
	1472	cccc is weighted by total ice mass (LESS GOOD than total ice surface area)
	1473	c rave(ig) = rave(ig) +
	1474	c & zq(ig,l,igcm_h2o_ice) *
	1475	c & (pplev(ig,l) - pplev(ig,l+1)) / g *
	1476	c & rice(ig,l) * (1.+nuice_ref)
[42]	1477	c Computing abs optical depth at 825 cm-1 in each
	1478	c layer to simulate NEW TES retrieval
	1479	Qabsice = min(
	1480	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1481	& )
	1482	opTES(ig,l)= 0.75 * Qabsice *
	1483	& zq(ig,l,igcm_h2o_ice) *
	1484	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1485	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1486	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1487	enddo
[669]	1488	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
	1489	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
[42]	1490	enddo
[669]	1491	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	1492	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
[42]	1493
[669]	1494	if (scavenging) then
	1495	Nccntot(:)= 0
	1496	Mccntot(:)= 0
	1497	rave(:)=0
	1498	do ig=1,ngrid
	1499	do l=1,nlayermx
	1500	Nccntot(ig) = Nccntot(ig) +
	1501	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
	1502	& *(pplev(ig,l) - pplev(ig,l+1)) / g
	1503	Mccntot(ig) = Mccntot(ig) +
	1504	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
	1505	& *(pplev(ig,l) - pplev(ig,l+1)) / g
	1506	cccc Column integrated effective ice radius
	1507	cccc is weighted by total ice surface area (BETTER than total ice mass)
	1508	rave(ig) = rave(ig) +
	1509	& tauscaling(ig) *
	1510	& zq(ig,l,igcm_ccn_number) *
	1511	& (pplev(ig,l) - pplev(ig,l+1)) / g *
	1512	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
	1513	enddo
	1514	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
	1515	& /max(pi*rave(ig),1.e-30) ! surface weight
	1516	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1517	enddo
	1518	endif ! of if (scavenging)
	1519
[42]	1520	endif ! of if (water)
	1521	endif ! of if (tracer)
	1522
	1523	c -----------------------------------------------------------------
	1524	c WSTATS: Saving statistics
	1525	c -----------------------------------------------------------------
	1526	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1527	c which can later be used to make the statistic files of the run:
	1528	c "stats") only possible in 3D runs !
	1529
[695]	1530	IF (callstats) THEN
[42]	1531
[695]	1532	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1533	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1534	call wstats(ngrid,"co2ice","CO2 ice cover",
[575]	1535	& "kg.m-2",2,co2ice)
[695]	1536	call wstats(ngrid,"fluxsurf_lw",
[575]	1537	& "Thermal IR radiative flux to surface","W.m-2",2,
	1538	& fluxsurf_lw)
[695]	1539	call wstats(ngrid,"fluxsurf_sw",
[575]	1540	& "Solar radiative flux to surface","W.m-2",2,
	1541	& fluxsurf_sw_tot)
[695]	1542	call wstats(ngrid,"fluxtop_lw",
[575]	1543	& "Thermal IR radiative flux to space","W.m-2",2,
	1544	& fluxtop_lw)
[695]	1545	call wstats(ngrid,"fluxtop_sw",
[575]	1546	& "Solar radiative flux to space","W.m-2",2,
	1547	& fluxtop_sw_tot)
[695]	1548	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1549	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1550	call wstats(ngrid,"v","Meridional (North-South) wind",
[575]	1551	& "m.s-1",3,zv)
[705]	1552	call wstats(ngrid,"w","Vertical (down-up) wind",
	1553	& "m.s-1",3,pw)
[695]	1554	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
[705]	1555	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
[226]	1556	c call wstats(ngrid,"q2",
	1557	c & "Boundary layer eddy kinetic energy",
	1558	c & "m2.s-2",3,q2)
	1559	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1560	c & emis)
	1561	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1562	c & 2,zstress)
	1563	c call wstats(ngrid,"sw_htrt","sw heat.rate",
	1564	c & "W.m-2",3,zdtsw)
	1565	c call wstats(ngrid,"lw_htrt","lw heat.rate",
	1566	c & "W.m-2",3,zdtlw)
[42]	1567
[226]	1568	if (tracer) then
	1569	if (water) then
[520]	1570	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1571	& *mugaz/mmol(igcm_h2o_vap)
	1572	call wstats(ngrid,"vmr_h2ovapor",
	1573	& "H2O vapor volume mixing ratio","mol/mol",
	1574	& 3,vmr)
	1575	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1576	& *mugaz/mmol(igcm_h2o_ice)
	1577	call wstats(ngrid,"vmr_h2oice",
	1578	& "H2O ice volume mixing ratio","mol/mol",
	1579	& 3,vmr)
[669]	1580	vmr=zqsat(1:ngridmx,1:nlayermx)
	1581	& *mugaz/mmol(igcm_h2o_vap)
	1582	call wstats(ngrid,"vmr_h2osat",
	1583	& "saturation volume mixing ratio","mol/mol",
	1584	& 3,vmr)
[226]	1585	call wstats(ngrid,"h2o_ice_s",
	1586	& "surface h2o_ice","kg/m2",
	1587	& 2,qsurf(1,igcm_h2o_ice))
[520]	1588	call wstats(ngrid,'albedo',
	1589	& 'albedo',
	1590	& '',2,albedo(1:ngridmx,1))
[226]	1591	call wstats(ngrid,"mtot",
	1592	& "total mass of water vapor","kg/m2",
	1593	& 2,mtot)
	1594	call wstats(ngrid,"icetot",
	1595	& "total mass of water ice","kg/m2",
	1596	& 2,icetot)
[520]	1597	call wstats(ngrid,"reffice",
	1598	& "Mean reff","m",
	1599	& 2,rave)
[669]	1600	call wstats(ngrid,"Nccntot",
[520]	1601	& "condensation nuclei","Nbr/m2",
[669]	1602	& 2,Nccntot)
	1603	call wstats(ngrid,"Mccntot",
	1604	& "condensation nuclei mass","kg/m2",
	1605	& 2,Mccntot)
[520]	1606	call wstats(ngrid,"rice",
	1607	& "Ice particle size","m",
	1608	& 3,rice)
[226]	1609	if (.not.activice) then
	1610	call wstats(ngrid,"tauTESap",
	1611	& "tau abs 825 cm-1","",
	1612	& 2,tauTES)
[520]	1613	else
	1614	call wstats(ngridmx,'tauTES',
	1615	& 'tau abs 825 cm-1',
	1616	& '',2,taucloudtes)
[226]	1617	endif
	1618
	1619	endif ! of if (water)
	1620
	1621	if (thermochem.or.photochem) then
	1622	do iq=1,nq
[576]	1623	if (noms(iq) .ne. "dust_mass" .and.
	1624	$ noms(iq) .ne. "dust_number" .and.
	1625	$ noms(iq) .ne. "ccn_mass" .and.
	1626	$ noms(iq) .ne. "ccn_number") then
[705]	1627	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	1628	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
	1629	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	1630	& rho(1:ngrid,1:nlayer)n_avog/
	1631	& (1000*mmol(iq))
[576]	1632	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1633	$ "Volume mixing ratio","mol/mol",3,vmr)
[705]	1634	! call wstats(ngrid,"rho_"//trim(noms(iq)),
	1635	! $ "Number density","cm-3",3,rhopart)
	1636	! call writediagfi(ngrid,"rho_"//trim(noms(iq)),
	1637	! $ "Number density","cm-3",3,rhopart)
[576]	1638	if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
	1639	$ (noms(iq).eq."o3")) then
	1640	call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
	1641	$ "Volume mixing ratio","mol/mol",3,vmr)
	1642	end if
	1643	do ig = 1,ngrid
	1644	colden(ig,iq) = 0.
	1645	end do
	1646	do l=1,nlayer
	1647	do ig=1,ngrid
	1648	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
	1649	$ *(pplev(ig,l)-pplev(ig,l+1))
	1650	$ 6.022e22/(mmol(iq)g)
	1651	end do
	1652	end do
	1653	call wstats(ngrid,"c_"//trim(noms(iq)),
	1654	$ "column","mol cm-2",2,colden(1,iq))
	1655	call writediagfi(ngrid,"c_"//trim(noms(iq)),
	1656	$ "column","mol cm-2",2,colden(1,iq))
	1657	end if ! of if (noms(iq) .ne. "dust_mass" ...)
	1658	end do ! of do iq=1,nq
[334]	1659	end if ! of if (thermochem.or.photochem)
[226]	1660
[334]	1661	end if ! of if (tracer)
[226]	1662
	1663	IF(lastcall) THEN
	1664	write (,) "Writing stats..."
	1665	call mkstats(ierr)
	1666	ENDIF
	1667
[42]	1668	ENDIF !if callstats
	1669
	1670	c (Store EOF for Mars Climate database software)
	1671	IF (calleofdump) THEN
	1672	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1673	ENDIF
	1674
[234]	1675
[226]	1676	#ifdef MESOSCALE
[234]	1677	!!!
	1678	!!! OUTPUT FIELDS
	1679	!!!
[81]	1680	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
[668]	1681	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
	1682	TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref)
[667]	1683	IF (tracer) THEN
[82]	1684	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
[608]	1685	icetot(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
[308]	1686	!! JF
	1687	IF (igcm_dust_mass .ne. 0) THEN
	1688	qsurfice_dust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
	1689	ENDIF
[81]	1690	IF (igcm_h2o_ice .ne. 0) THEN
	1691	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
	1692	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	1693	. * mugaz / mmol(igcm_h2o_ice)
	1694	ENDIF
[308]	1695	!! Dust quantity integration along the vertical axe
	1696	dustot(:)=0
[698]	1697	IF (igcm_dust_mass .ne. 0) THEN
[308]	1698	do ig=1,ngrid
	1699	do l=1,nlayermx
	1700	dustot(ig) = dustot(ig) +
	1701	& zq(ig,l,igcm_dust_mass)
	1702	& * (pplev(ig,l) - pplev(ig,l+1)) / g
	1703	enddo
	1704	enddo
[667]	1705	ENDIF
[698]	1706	ENDIF
[550]	1707	!! TAU water ice as seen by TES
	1708	if (activice) tauTES = taucloudtes
[226]	1709	c AUTOMATICALLY GENERATED FROM REGISTRY
[42]	1710	#include "fill_save.inc"
[226]	1711	#else
[528]	1712	#ifndef MESOINI
[42]	1713
[226]	1714	c ==========================================================
	1715	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1716	c any variable for diagnostic (output with period
	1717	c "ecritphy", set in "run.def")
	1718	c ==========================================================
	1719	c WRITEDIAGFI can ALSO be called from any other subroutines
	1720	c for any variables !!
	1721	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1722	c & emis)
[520]	1723	c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
	1724	c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
[226]	1725	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1726	& tsurf)
	1727	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1728	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1729	& co2ice)
[284]	1730
[226]	1731	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1732	c & "K",2,zt(1,7))
[284]	1733	c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1734	c & fluxsurf_lw)
	1735	c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1736	c & fluxsurf_sw_tot)
	1737	c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1738	c & fluxtop_lw)
	1739	c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1740	c & fluxtop_sw_tot)
[520]	1741	c call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1742	c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1743	c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1744	c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1745	c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
[226]	1746	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
[520]	1747	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
[226]	1748	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
	1749	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1750	c & zstress)
	1751	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	1752	c & 'w.m-2',3,zdtsw)
	1753	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	1754	c & 'w.m-2',3,zdtlw)
[520]	1755	if (.not.activice) then
	1756	CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1757	& 'tau abs 825 cm-1',
	1758	& '',2,tauTES)
	1759	else
	1760	CALL WRITEDIAGFI(ngridmx,'tauTES',
	1761	& 'tau abs 825 cm-1',
	1762	& '',2,taucloudtes)
	1763	endif
[42]	1764
[528]	1765	#else
	1766	!!! this is to ensure correct initialisation of mesoscale model
	1767	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1768	& tsurf)
	1769	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1770	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1771	& co2ice)
	1772	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1773	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1774	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
[299]	1775	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1776	& emis)
	1777	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
	1778	& "K",3,tsoil)
	1779	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
	1780	& "K",3,inertiedat)
	1781	#endif
[42]	1782
[299]	1783
[226]	1784	c ----------------------------------------------------------
	1785	c Outputs of the CO2 cycle
	1786	c ----------------------------------------------------------
[42]	1787
[226]	1788	if (tracer.and.(igcm_co2.ne.0)) then
	1789	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	1790	! & "kg/kg",2,zq(1,1,igcm_co2))
	1791	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	1792	! & "kg/kg",3,zq(1,1,igcm_co2))
	1793
	1794	! Compute co2 column
[286]	1795	co2col(:)=0
[226]	1796	do l=1,nlayermx
	1797	do ig=1,ngrid
	1798	co2col(ig)=co2col(ig)+
	1799	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
	1800	enddo
	1801	enddo
	1802	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	1803	& co2col)
	1804	endif ! of if (tracer.and.(igcm_co2.ne.0))
[42]	1805
[226]	1806	c ----------------------------------------------------------
	1807	c Outputs of the water cycle
	1808	c ----------------------------------------------------------
	1809	if (tracer) then
	1810	if (water) then
	1811
[299]	1812	#ifdef MESOINI
	1813	!!!! waterice = q01, voir readmeteo.F90
	1814	call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
	1815	& 'kg/kg',3,
	1816	& zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
	1817	!!!! watervapor = q02, voir readmeteo.F90
	1818	call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap),
	1819	& 'kg/kg',3,
	1820	& zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
	1821	!!!! surface waterice qsurf02 (voir readmeteo)
	1822	call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer',
	1823	& 'kg.m-2',2,
	1824	& qsurf(1:ngridmx,igcm_h2o_ice))
	1825	#endif
	1826
[226]	1827	CALL WRITEDIAGFI(ngridmx,'mtot',
	1828	& 'total mass of water vapor',
	1829	& 'kg/m2',2,mtot)
	1830	CALL WRITEDIAGFI(ngridmx,'icetot',
	1831	& 'total mass of water ice',
	1832	& 'kg/m2',2,icetot)
[411]	1833	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1834	c & *mugaz/mmol(igcm_h2o_ice)
	1835	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
	1836	c & 'mol/mol',3,vmr)
	1837	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1838	c & *mugaz/mmol(igcm_h2o_vap)
	1839	c call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
	1840	c & 'mol/mol',3,vmr)
[669]	1841	CALL WRITEDIAGFI(ngridmx,'reffice',
	1842	& 'Mean reff',
	1843	& 'm',2,rave)
	1844	CALL WRITEDIAGFI(ngrid,"Nccntot",
	1845	& "condensation nuclei","Nbr/m2",
	1846	& 2,Nccntot)
	1847	c CALL WRITEDIAGFI(ngrid,"Mccntot",
	1848	c & "mass condensation nuclei","kg/m2",
	1849	c & 2,Mccntot)
[411]	1850	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
	1851	c & 'm',3,rice)
[226]	1852	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
	1853	& 'surface h2o_ice',
	1854	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
[358]	1855	c CALL WRITEDIAGFI(ngridmx,'albedo',
	1856	c & 'albedo',
	1857	c & '',2,albedo(1:ngridmx,1))
[226]	1858	endif !(water)
	1859
	1860
	1861	if (water.and..not.photochem) then
	1862	iq=nq
	1863	c write(str2(1:2),'(i2.2)') iq
	1864	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
	1865	c & 'kg.m-2',2,zdqscloud(1,iq))
	1866	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
	1867	c & 'kg/kg',3,zdqchim(1,1,iq))
	1868	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
	1869	c & 'kg/kg',3,zdqdif(1,1,iq))
	1870	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
	1871	c & 'kg/kg',3,zdqadj(1,1,iq))
	1872	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
	1873	c & 'kg/kg',3,zdqc(1,1,iq))
	1874	endif !(water.and..not.photochem)
	1875	endif
	1876
	1877	c ----------------------------------------------------------
	1878	c Outputs of the dust cycle
	1879	c ----------------------------------------------------------
	1880
	1881	c call WRITEDIAGFI(ngridmx,'tauref',
	1882	c & 'Dust ref opt depth','NU',2,tauref)
	1883
	1884	if (tracer.and.(dustbin.ne.0)) then
	1885	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
	1886	if (doubleq) then
[411]	1887	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
	1888	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
	1889	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
	1890	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
[226]	1891	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
	1892	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
[411]	1893	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
	1894	c & 'kg.m-2.s-1',2,zdqssed(1,1))
	1895	c call WRITEDIAGFI(ngridmx,'dqsdif','diffusion',
	1896	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
	1897	c call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
	1898	c & 'm',3,rdust*ref_r0)
	1899	c call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	1900	c & 'kg/kg',3,pq(1,1,igcm_dust_mass))
	1901	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
	1902	c & 'part/kg',3,pq(1,1,igcm_dust_number))
[299]	1903	#ifdef MESOINI
[550]	1904	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	1905	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
[299]	1906	call WRITEDIAGFI(ngridmx,'dustN','Dust number',
[520]	1907	& 'part/kg',3,pq(1,1,igcm_dust_number))
[550]	1908	call WRITEDIAGFI(ngridmx,'ccn','Nuclei mass mr',
	1909	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
	1910	call WRITEDIAGFI(ngridmx,'ccnN','Nuclei number',
	1911	& 'part/kg',3,pq(1,1,igcm_ccn_number))
[299]	1912	#endif
[226]	1913	else
	1914	do iq=1,dustbin
	1915	write(str2(1:2),'(i2.2)') iq
	1916	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
	1917	& 'kg/kg',3,zq(1,1,iq))
	1918	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
	1919	& 'kg.m-2',2,qsurf(1,iq))
	1920	end do
	1921	endif ! (doubleq)
[358]	1922
	1923	if (scavenging) then
[520]	1924	c call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
	1925	c & 'kg/kg',3,pq(1,1,igcm_ccn_mass))
	1926	c call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
	1927	c & 'part/kg',3,pq(1,1,igcm_ccn_number))
[358]	1928	endif ! (scavenging)
	1929
[226]	1930	c if (submicron) then
	1931	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
	1932	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
	1933	c endif ! (submicron)
	1934	end if ! (tracer.and.(dustbin.ne.0))
	1935
[705]	1936
[226]	1937	c ----------------------------------------------------------
[705]	1938	c Thermospheric outputs
[267]	1939	c ----------------------------------------------------------
[705]	1940
	1941	if(callthermos) then
	1942
	1943	call WRITEDIAGFI(ngridmx,"q15um","15 um cooling","K/s",
	1944	$ 3,zdtnlte)
	1945	call WRITEDIAGFI(ngridmx,"quv","UV heating","K/s",
	1946	$ 3,zdteuv)
	1947	call WRITEDIAGFI(ngridmx,"cond","Thermal conduction","K/s",
	1948	$ 3,zdtconduc)
	1949	call WRITEDIAGFI(ngridmx,"qnir","NIR heating","K/s",
	1950	$ 3,zdtnirco2)
	1951
	1952	endif !(callthermos)
	1953
	1954	c ----------------------------------------------------------
	1955	c ----------------------------------------------------------
[267]	1956	c PBL OUTPUS
	1957	c ----------------------------------------------------------
	1958	c ----------------------------------------------------------
	1959
	1960	c ----------------------------------------------------------
[226]	1961	c Outputs of thermals
	1962	c ----------------------------------------------------------
	1963	if (calltherm) then
	1964
	1965	! call WRITEDIAGFI(ngrid,'dtke',
	1966	! & 'tendance tke thermiques','m2/s2',
	1967	! & 3,dtke_th)
	1968	! call WRITEDIAGFI(ngrid,'d_u_ajs',
	1969	! & 'tendance u thermiques','m/s',
	1970	! & 3,pdu_th*ptimestep)
	1971	! call WRITEDIAGFI(ngrid,'d_v_ajs',
	1972	! & 'tendance v thermiques','m/s',
	1973	! & 3,pdv_th*ptimestep)
	1974	! if (tracer) then
	1975	! if (nq .eq. 2) then
	1976	! call WRITEDIAGFI(ngrid,'deltaq_th',
	1977	! & 'delta q thermiques','kg/kg',
	1978	! & 3,ptimestep*pdq_th(:,:,2))
	1979	! endif
	1980	! endif
	1981
[523]	1982	call WRITEDIAGFI(ngridmx,'zmax_th',
	1983	& 'hauteur du thermique','m',
	1984	& 2,zmax_th)
[226]	1985	call WRITEDIAGFI(ngridmx,'hfmax_th',
	1986	& 'maximum TH heat flux','K.m/s',
	1987	& 2,hfmax_th)
[528]	1988	call WRITEDIAGFI(ngridmx,'wstar',
[226]	1989	& 'maximum TH vertical velocity','m/s',
[528]	1990	& 2,wstar)
[226]	1991
	1992	endif
	1993
	1994	c ----------------------------------------------------------
[267]	1995	c ----------------------------------------------------------
	1996	c END OF PBL OUTPUS
	1997	c ----------------------------------------------------------
	1998	c ----------------------------------------------------------
	1999
	2000
	2001	c ----------------------------------------------------------
[226]	2002	c Output in netcdf file "diagsoil.nc" for subterranean
	2003	c variables (output every "ecritphy", as for writediagfi)
	2004	c ----------------------------------------------------------
	2005
	2006	! Write soil temperature
	2007	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	2008	! & 3,tsoil)
	2009	! Write surface temperature
	2010	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	2011	! & 2,tsurf)
	2012
	2013	c ==========================================================
	2014	c END OF WRITEDIAGFI
	2015	c ==========================================================
	2016	#endif
	2017
[42]	2018	ELSE ! if(ngrid.eq.1)
	2019
[627]	2020	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
	2021	& zls*180./pi,odpref,tauref
[42]	2022	c ----------------------------------------------------------------------
	2023	c Output in grads file "g1d" (ONLY when using testphys1d)
	2024	c (output at every X physical timestep)
	2025	c ----------------------------------------------------------------------
	2026	c
	2027	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	2028	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	2029	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	2030
	2031	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	2032	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	2033	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	2034	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	2035
[226]	2036	! THERMALS STUFF 1D
	2037	if(calltherm) then
	2038
	2039	call WRITEDIAGFI(ngridmx,'lmax_th',
	2040	& 'hauteur du thermique','point',
	2041	& 0,lmax_th_out)
[528]	2042	call WRITEDIAGFI(ngridmx,'zmax_th',
	2043	& 'hauteur du thermique','m',
	2044	& 0,zmax_th)
[226]	2045	call WRITEDIAGFI(ngridmx,'hfmax_th',
	2046	& 'maximum TH heat flux','K.m/s',
	2047	& 0,hfmax_th)
[528]	2048	call WRITEDIAGFI(ngridmx,'wstar',
[226]	2049	& 'maximum TH vertical velocity','m/s',
[528]	2050	& 0,wstar)
[226]	2051
	2052	co2col(:)=0.
	2053	if (tracer) then
	2054	do l=1,nlayermx
	2055	do ig=1,ngrid
	2056	co2col(ig)=co2col(ig)+
	2057	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
	2058	enddo
	2059	enddo
	2060
	2061	end if
	2062	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
	2063	& ,'kg/m-2',0,co2col)
	2064	endif
	2065	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
	2066	& ,'m/s',1,pw)
	2067	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
	2068	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
	2069	& tsurf)
[277]	2070	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
	2071	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
[226]	2072
	2073	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
	2074	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
[358]	2075	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
[544]	2076	! call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
	2077	! & "K.s-1",1,dtrad/zpopsk)
	2078	! call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	2079	! & 'w.m-2',1,zdtsw/zpopsk)
	2080	! call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	2081	! & 'w.m-2',1,zdtlw/zpopsk)
	2082
[226]	2083	! or output in diagfi.nc (for testphys1d)
	2084	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
	2085	call WRITEDIAGFI(ngridmx,'temp','Temperature',
	2086	& 'K',1,zt)
[358]	2087
[226]	2088	if(tracer) then
	2089	c CALL writeg1d(ngrid,1,tau,'tau','SI')
	2090	do iq=1,nq
	2091	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	2092	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
	2093	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	2094	end do
[358]	2095	if (doubleq) then
	2096	call WRITEDIAGFI(ngridmx,'rdust','rdust',
	2097	& 'm',1,rdust)
	2098	endif
[226]	2099	end if
[358]	2100
[520]	2101	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
[358]	2102	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[520]	2103	IF (water) THEN
	2104	CALL WRITEDIAGFI(ngridmx,'tauTESap',
[358]	2105	& 'tau abs 825 cm-1',
[520]	2106	& '',0,tauTES)
[358]	2107
[520]	2108	CALL WRITEDIAGFI(ngridmx,'tauTES',
	2109	& 'tau abs 825 cm-1',
	2110	& '',0,taucloudtes)
	2111
	2112	mtot = 0
	2113	icetot = 0
	2114	h2otot = qsurf(1,igcm_h2o_ice)
	2115	rave = 0
[358]	2116	do l=1,nlayer
[520]	2117	mtot = mtot + zq(1,l,igcm_h2o_vap)
	2118	& * (pplev(1,l) - pplev(1,l+1)) / g
	2119	icetot = icetot + zq(1,l,igcm_h2o_ice)
	2120	& * (pplev(1,l) - pplev(1,l+1)) / g
[669]	2121	cccc Column integrated effective ice radius
	2122	cccc is weighted by total ice surface area (BETTER)
	2123	rave = rave + tauscaling(ig) *
	2124	& zq(1,l,igcm_ccn_number) *
	2125	& (pplev(1,l) - pplev(1,l+1)) / g *
	2126	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
	2127	cccc Column integrated effective ice radius
	2128	cccc is weighted by total ice mass (LESS GOOD)
	2129	c rave = rave + zq(1,l,igcm_h2o_ice)
	2130	c & * (pplev(1,l) - pplev(1,l+1)) / g
	2131	c & * rice(1,l) * (1.+nuice_ref)
[358]	2132	end do
[669]	2133	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
	2134	c rave=rave/max(icetot,1.e-30) ! mass weight
[520]	2135	h2otot = h2otot+mtot+icetot
	2136
	2137
	2138	if (scavenging) then
[669]	2139	Nccntot= 0
[520]	2140	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	2141	do l=1,nlayermx
[669]	2142	Nccntot = Nccntot +
[520]	2143	& zq(1,l,igcm_ccn_number)*tauscaling(1)
	2144	& *(pplev(1,l) - pplev(1,l+1)) / g
	2145	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
	2146	satu(1,l) = (max(satu(1,l),float(1))-1)
	2147	! & * zq(1,l,igcm_h2o_vap) *
	2148	! & (pplev(1,l) - pplev(1,l+1)) / g
	2149	enddo
	2150
[669]	2151	CALL WRITEDIAGFI(ngridmx,'Nccntot',
	2152	& 'Nccntot',
	2153	& 'nbr/m2',0,Nccntot)
[520]	2154	endif
	2155
	2156
	2157	CALL WRITEDIAGFI(ngridmx,'h2otot',
	2158	& 'h2otot',
	2159	& 'kg/m2',0,h2otot)
	2160	CALL WRITEDIAGFI(ngridmx,'mtot',
	2161	& 'mtot',
	2162	& 'kg/m2',0,mtot)
	2163	CALL WRITEDIAGFI(ngridmx,'icetot',
	2164	& 'icetot',
	2165	& 'kg/m2',0,icetot)
	2166	CALL WRITEDIAGFI(ngridmx,'reffice',
	2167	& 'reffice',
	2168	& 'm',0,rave)
[411]	2169
	2170
	2171	do iq=1,nq
	2172	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_s',
	2173	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
	2174	end do
	2175
[358]	2176
[520]	2177	call WRITEDIAGFI(ngridmx,'zdqsed_dustq','sedimentation q',
	2178	& 'kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
	2179	call WRITEDIAGFI(ngridmx,'zdqsed_dustN','sedimentation N',
	2180	& 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_dust_number))
[358]	2181
[633]	2182	call WRITEDIAGFI(ngridmx,'zdqcloud_ice','cloud ice',
	2183	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
	2184	call WRITEDIAGFI(ngridmx,'zdqcloud_vap','cloud vap',
	2185	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
	2186	call WRITEDIAGFI(ngridmx,'zdqcloud','cloud ice',
	2187	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
	2188	& +zdqcloud(1,:,igcm_h2o_vap))
[411]	2189
[633]	2190
[520]	2191	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
[411]	2192	& rice)
[520]	2193	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
	2194	& satu)
	2195	ENDIF ! of IF (water)
[411]	2196
[358]	2197	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2198	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2199
	2200
[226]	2201	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
	2202
	2203	do l=2,nlayer-1
	2204	tmean=zt(1,l)
	2205	if(zt(1,l).ne.zt(1,l-1))
	2206	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	2207	zlocal(l)= zlocal(l-1)
	2208	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
	2209	enddo
	2210	zlocal(nlayer)= zlocal(nlayer-1)-
	2211	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
	2212	& rnew(1,nlayer)*tmean/g
	2213
[42]	2214	END IF ! if(ngrid.ne.1)
	2215
	2216	icount=icount+1
	2217	RETURN
	2218	END

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

Download in other formats:

Original Format