Context Navigation

physiq.F @ 735

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

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

Download in other formats:

Original Format