Context Navigation

physiq.F @ 887

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

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

Download in other formats:

Original Format