Context Navigation

physiq.F @ 551

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

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

Download in other formats:

Original Format