Context Navigation

physiq.F @ 758

Last change on this file since 758 was 758, checked in by emillour, 12 years ago

Mars GCM:

Some minor changes and adaptations for the MCDv5 runs:

removed the 'nice' instruction from run0 script
update dissipation factors in inidissip.F (fac_mid=3,fac_up=30,startalt=70)
add some outputs in physiq.F
put example def files used for MCDv5 simulations in deftank: callphys.def.MCD5 run.def.64x48x49.MCD5 traceur.def.MCD5

File size: 87.8 KB

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

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

Download in other formats:

Original Format