Context Navigation

physiq.F @ 185

Last change on this file since 185 was 185, checked in by acolaitis, 14 years ago

17/06/2011 == AC

Added new settings for the Martian thermals from new LES observations
Revamped thermcell's module variables to allow it's removal
Minor changes in physiq and meso_physiq for the call to thermals
Switched from dynamic to static memory allocation for all thermals variable

to gain computation speed

Thermals now output maximum speed, maximum heat flux, and maximum height in thermal plumes

File size: 57.5 KB

Rev	Line
[38]	1	SUBROUTINE physiq(ngrid,nlayer,nq,
	2	$ firstcall,lastcall,
	3	$ pday,ptime,ptimestep,
	4	$ pplev,pplay,pphi,
	5	$ pu,pv,pt,pq,
	6	$ pw,
	7	$ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
	8
	9	IMPLICIT NONE
	10	c=======================================================================
	11	c
	12	c subject:
	13	c --------
	14	c
	15	c Organisation of the physical parametrisations of the LMD
	16	c martian atmospheric general circulation model.
	17	c
	18	c The GCM can be run without or with tracer transport
	19	c depending on the value of Logical "tracer" in file "callphys.def"
	20	c Tracers may be water vapor, ice OR chemical species OR dust particles
	21	c
	22	c SEE comments in initracer.F about numbering of tracer species...
	23	c
	24	c It includes:
	25	c
	26	c 1. Initialization:
	27	c 1.1 First call initializations
	28	c 1.2 Initialization for every call to physiq
	29	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	30	c 2. Compute radiative transfer tendencies
	31	c (longwave and shortwave) for CO2 and aerosols.
	32	c 3. Gravity wave and subgrid scale topography drag :
	33	c 4. Vertical diffusion (turbulent mixing):
	34	c 5. Convective adjustment
	35	c 6. Condensation and sublimation of carbon dioxide.
	36	c 7. TRACERS :
	37	c 7a. water and water ice
	38	c 7b. call for photochemistry when tracers are chemical species
	39	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
	40	c 7d. updates (CO2 pressure variations, surface budget)
	41	c 8. Contribution to tendencies due to thermosphere
	42	c 9. Surface and sub-surface temperature calculations
	43	c 10. Write outputs :
	44	c - "startfi", "histfi" (if it's time)
	45	c - Saving statistics (if "callstats = .true.")
	46	c - Dumping eof (if "calleofdump = .true.")
	47	c - Output any needed variables in "diagfi"
	48	c 11. Diagnostic: mass conservation of tracers
	49	c
	50	c author:
	51	c -------
	52	c Frederic Hourdin 15/10/93
	53	c Francois Forget 1994
	54	c Christophe Hourdin 02/1997
	55	c Subroutine completly rewritten by F.Forget (01/2000)
	56	c Introduction of the photochemical module: S. Lebonnois (11/2002)
	57	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
	58	c Water ice clouds: Franck Montmessin (update 06/2003)
	59	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	60	c Nb: See callradite.F for more information.
	61	c
	62	c arguments:
	63	c ----------
	64	c
	65	c input:
	66	c ------
	67	c ecri period (in dynamical timestep) to write output
	68	c ngrid Size of the horizontal grid.
	69	c All internal loops are performed on that grid.
	70	c nlayer Number of vertical layers.
	71	c nq Number of advected fields
	72	c firstcall True at the first call
	73	c lastcall True at the last call
	74	c pday Number of days counted from the North. Spring
	75	c equinoxe.
	76	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	77	c ptimestep timestep (s)
	78	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	79	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	80	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	81	c pu(ngrid,nlayer) u component of the wind (ms-1)
	82	c pv(ngrid,nlayer) v component of the wind (ms-1)
	83	c pt(ngrid,nlayer) Temperature (K)
	84	c pq(ngrid,nlayer,nq) Advected fields
	85	c pudyn(ngrid,nlayer) \
	86	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	87	c ptdyn(ngrid,nlayer) / corresponding variables
	88	c pqdyn(ngrid,nlayer,nq) /
	89	c pw(ngrid,?) vertical velocity
	90	c
	91	c output:
	92	c -------
	93	c
	94	c pdu(ngrid,nlayermx) \
	95	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
	96	c pdt(ngrid,nlayermx) / variables due to physical processes.
	97	c pdq(ngrid,nlayermx,nqmx) /
	98	c pdpsrf(ngrid) /
	99	c tracerdyn call tracer in dynamical part of GCM ?
	100
	101	c
	102	c=======================================================================
	103	c
	104	c 0. Declarations :
	105	c ------------------
	106
	107	#include "dimensions.h"
	108	#include "dimphys.h"
	109	#include "comgeomfi.h"
	110	#include "surfdat.h"
	111	#include "comsoil.h"
	112	#include "comdiurn.h"
	113	#include "callkeys.h"
	114	#include "comcstfi.h"
	115	#include "planete.h"
	116	#include "comsaison.h"
	117	#include "control.h"
	118	#include "dimradmars.h"
	119	#include "comg1d.h"
	120	#include "tracer.h"
	121	#include "nlteparams.h"
	122
	123	#include "chimiedata.h"
	124	#include "watercap.h"
	125	#include "param.h"
	126	#include "param_v3.h"
	127	#include "conc.h"
	128
	129	#include "netcdf.inc"
	130
	131
	132
	133	c Arguments :
	134	c -----------
	135
	136	c inputs:
	137	c -------
	138	INTEGER ngrid,nlayer,nq
	139	REAL ptimestep
	140	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	141	REAL pphi(ngridmx,nlayer)
	142	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	143	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	144	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
	145	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	146	LOGICAL firstcall,lastcall
	147
	148	REAL pday
	149	REAL ptime
	150	logical tracerdyn
	151
	152	c outputs:
	153	c --------
	154	c physical tendencies
	155	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	156	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	157	REAL pdpsrf(ngridmx) ! surface pressure tendency
	158
	159
	160	c Local saved variables:
	161	c ----------------------
	162	c aerosol (dust or ice) extinction optical depth at reference wavelength
	163	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
	164	c aerosol optical properties :
	165	REAL aerosol(ngridmx,nlayermx,naerkind)
	166
	167	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	168	INTEGER icount ! counter of calls to physiq during the run.
	169	REAL tsurf(ngridmx) ! Surface temperature (K)
	170	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	171	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
	172	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
	173	REAL emis(ngridmx) ! Thermal IR surface emissivity
	174	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	175	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	176	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	177	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	178	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
	179	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	180	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	181	INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic)
	182
	183	c Variables used by the water ice microphysical scheme:
	184	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
	185	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
	186	! of the size distribution
	187	c Albedo of deposited surface ice
[82]	188	!!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
	189	REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
[38]	190
	191	SAVE day_ini, icount
	192	SAVE aerosol, tsurf,tsoil
	193	SAVE co2ice,albedo,emis, q2
	194	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
	195	SAVE ig_vl1
	196
	197	REAL stephan
	198	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	199	SAVE stephan
	200
	201	c Local variables :
	202	c -----------------
	203
	204	REAL CBRT
	205	EXTERNAL CBRT
	206
	207	CHARACTER*80 fichier
[185]	208	INTEGER l,ig,ierr,igout,iq,i, tapphys,k
	209	LOGICAL end_of_file
[38]	210
	211	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
	212	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
	213	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
	214	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
	215	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
	216	! (used if active=F)
	217	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	218	REAL zls ! solar longitude (rad)
	219	REAL zday ! date (time since Ls=0, in martian days)
	220	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	221	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
	222	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
	223
	224	c Tendancies due to various processes:
	225	REAL dqsurf(ngridmx,nqmx)
	226	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	227	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
	228	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
	229	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
	230	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
	231	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
	232	REAL zdtsurf(ngridmx) ! (K/s)
	233	REAL zdtcloud(ngridmx,nlayermx)
	234	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	235	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	236	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	237	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	238	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	239	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	240	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	241	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	242
	243	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
	244	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
	245	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
	246	REAL zdqadj(ngridmx,nlayermx,nqmx)
	247	REAL zdqc(ngridmx,nlayermx,nqmx)
	248	REAL zdqcloud(ngridmx,nlayermx,nqmx)
	249	REAL zdqscloud(ngridmx,nqmx)
	250	REAL zdqchim(ngridmx,nlayermx,nqmx)
	251	REAL zdqschim(ngridmx,nqmx)
	252
	253	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	254	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	255	REAL zdumolvis(ngridmx,nlayermx)
	256	REAL zdvmolvis(ngridmx,nlayermx)
	257	real zdqmoldiff(ngridmx,nlayermx,nqmx)
	258
	259	c Local variable for local intermediate calcul:
	260	REAL zflubid(ngridmx)
	261	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	262	REAL zdum1(ngridmx,nlayermx)
	263	REAL zdum2(ngridmx,nlayermx)
	264	REAL ztim1,ztim2,ztim3, z1,z2
	265	REAL ztime_fin
	266	REAL zdh(ngridmx,nlayermx)
	267	INTEGER length
	268	PARAMETER (length=100)
	269
	270	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	271	c -----------------------------------------------------------------------------
	272	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	273	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	274	REAL zq(ngridmx,nlayermx,nqmx)
	275	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	276	character*2 str2
	277	character*5 str5
	278	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
	279	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
	280	! (particules kg-1)
[83]	281	SAVE ccn !! in case iradia != 1
[38]	282	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
	283	real qtot1,qtot2 ! total aerosol mass
	284	integer igmin, lmin
	285	logical tdiag
	286
	287	real co2col(ngridmx) ! CO2 column
	288	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
	289	REAL zstress(ngrid), cd
	290	real hco2(nqmx),tmean, zlocal(nlayermx)
	291	real rho(ngridmx,nlayermx) ! density
	292	real vmr(ngridmx,nlayermx) ! volume mixing ratio
	293	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
	294	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
	295	REAL rave(ngridmx) ! Mean water ice effective radius (m)
	296	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
	297	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
	298	REAL Qabsice ! Water ice absorption coefficient
	299
	300
	301	REAL time_phys
	302
[161]	303	c Variables from thermal
	304
[185]	305	REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
	306	REAL wmax_th(ngridmx)
	307	REAL ,SAVE :: hfmax_th(ngridmx)
	308	REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
	309	REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
	310	INTEGER lmax_th(ngridmx)
	311	REAL dtke_th(ngridmx,nlayermx+1)
	312	REAL dummycol(ngridmx)
[173]	313
[38]	314	c=======================================================================
	315
	316	c 1. Initialisation:
	317	c -----------------
	318
	319	c 1.1 Initialisation only at first call
	320	c ---------------------------------------
	321	IF (firstcall) THEN
	322
	323	c variables set to 0
	324	c ~~~~~~~~~~~~~~~~~~
	325	call zerophys(ngridnlayernaerkind,aerosol)
	326	call zerophys(ngrid*nlayer,dtrad)
	327	call zerophys(ngrid,fluxrad)
	328
	329	c read startfi
	330	c ~~~~~~~~~~~~
	331
	332	! Read netcdf initial physical parameters.
	333	CALL phyetat0 ("startfi.nc",0,0,
	334	& nsoilmx,nq,
	335	& day_ini,time_phys,
	336	& tsurf,tsoil,emis,q2,qsurf,co2ice)
	337
	338	if (pday.ne.day_ini) then
	339	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	340	& "physics and dynamics"
	341	write(,) "dynamics day: ",pday
	342	write(,) "physics day: ",day_ini
	343	stop
	344	endif
	345
	346	write (,) 'In physiq day_ini =', day_ini
	347
	348	c Initialize albedo and orbital calculation
	349	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	350	CALL surfini(ngrid,co2ice,qsurf,albedo)
	351	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	352
	353	c initialize soil
	354	c ~~~~~~~~~~~~~~~
	355	IF (callsoil) THEN
	356	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	357	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	358	ELSE
	359	PRINT*,
	360	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	361	DO ig=1,ngrid
	362	capcal(ig)=1.e5
	363	fluxgrd(ig)=0.
	364	ENDDO
	365	ENDIF
	366	icount=1
	367
	368
	369	c initialize tracers
	370	c ~~~~~~~~~~~~~~~~~~
	371	tracerdyn=tracer
	372	IF (tracer) THEN
	373	CALL initracer(qsurf,co2ice)
	374	ENDIF ! end tracer
	375
	376	c Determining gridpoint near Viking Lander 1 (used for diagnostic only)
	377	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	378
	379	if(ngrid.ne.1) then
	380	latvl1= 22.27
	381	lonvl1= -47.94
	382	ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim +
	383	& int(1.5+(lonvl1+180)*iim/360.)
	384	write(,) 'Viking Lander 1 GCM point: lat,lon',
	385	& lati(ig_vl1)180/pi, long(ig_vl1)180/pi
	386	end if
	387
	388	c Initialize thermospheric parameters
	389	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	390
	391	if (callthermos) call param_read
	392
	393	c Initialize R and Cp as constant
	394
	395	if (.not.callthermos .and. .not.photochem) then
	396	do l=1,nlayermx
	397	do ig=1,ngridmx
	398	rnew(ig,l)=r
	399	cpnew(ig,l)=cpp
	400	mmean(ig,l)=mugaz
	401	enddo
	402	enddo
	403	endif
	404
	405	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
	406	write(,)"physiq: water_param Surface ice alb:",alb_surfice
	407	ENDIF
	408
	409	ENDIF ! (end of "if firstcall")
	410
	411	c ---------------------------------------------------
	412	c 1.2 Initializations done at every physical timestep:
	413	c ---------------------------------------------------
	414	c
	415	IF (ngrid.NE.ngridmx) THEN
	416	PRINT*,'STOP in PHYSIQ'
	417	PRINT*,'Probleme de dimensions :'
	418	PRINT*,'ngrid = ',ngrid
	419	PRINT*,'ngridmx = ',ngridmx
	420	STOP
	421	ENDIF
	422
	423	c Initialize various variables
	424	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	425	call zerophys(ngrid*nlayer, pdv)
	426	call zerophys(ngrid*nlayer, pdu)
	427	call zerophys(ngrid*nlayer, pdt)
	428	call zerophys(ngridnlayernq, pdq)
	429	call zerophys(ngrid, pdpsrf)
	430	call zerophys(ngrid, zflubid)
	431	call zerophys(ngrid, zdtsurf)
	432	call zerophys(ngrid*nq, dqsurf)
	433	igout=ngrid/2+1
	434
	435
	436	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	437
	438	c Compute Solar Longitude (Ls) :
	439	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	440	if (season) then
	441	call solarlong(zday,zls)
	442	else
	443	call solarlong(float(day_ini),zls)
	444	end if
	445
	446	c Compute geopotential at interlayers
	447	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	448	c ponderation des altitudes au niveau des couches en dp/p
	449
	450	DO l=1,nlayer
	451	DO ig=1,ngrid
	452	zzlay(ig,l)=pphi(ig,l)/g
	453	ENDDO
	454	ENDDO
	455	DO ig=1,ngrid
	456	zzlev(ig,1)=0.
	457	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	458	ENDDO
	459	DO l=2,nlayer
	460	DO ig=1,ngrid
	461	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	462	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	463	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	464	ENDDO
	465	ENDDO
	466
	467
	468	! Potential temperature calculation not the same in physiq and dynamic
	469
	470	c Compute potential temperature
	471	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	472	DO l=1,nlayer
	473	DO ig=1,ngrid
	474	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	475	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	476	ENDDO
	477	ENDDO
	478
	479	c-----------------------------------------------------------------------
	480	c 1.2.5 Compute mean mass, cp, and R
	481	c --------------------------------
	482
	483	if(photochem.or.callthermos) then
	484	call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
	485	endif
	486	c-----------------------------------------------------------------------
	487	c 2. Compute radiative tendencies :
	488	c------------------------------------
	489
	490
	491	IF (callrad) THEN
	492	IF( MOD(icount-1,iradia).EQ.0) THEN
	493
	494	c Local Solar zenith angle
	495	c ~~~~~~~~~~~~~~~~~~~~~~~~
	496	CALL orbite(zls,dist_sol,declin)
	497
	498	IF(diurnal) THEN
	499	ztim1=SIN(declin)
	500	ztim2=COS(declin)COS(2.pi*(zday-.5))
	501	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	502
	503	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	504	s ztim1,ztim2,ztim3, mu0,fract)
	505
	506	ELSE
	507	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
	508	ENDIF
	509
	510	c NLTE cooling from CO2 emission
	511	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	512
	513	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
	514
	515	c Find number of layers for LTE radiation calculations
	516	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
	517	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
	518
	519	c Note: Dustopacity.F has been transferred to callradite.F
	520
	521	c Call main radiative transfer scheme
	522	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	523	c Transfer through CO2 (except NIR CO2 absorption)
	524	c and aerosols (dust and water ice)
	525
	526	c Radiative transfer
	527	c ------------------
[161]	528
[38]	529	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
	530	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
	531	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
	532	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
	533
	534	c CO2 near infrared absorption
	535	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	536	call zerophys(ngrid*nlayer,zdtnirco2)
	537	if (callnirco2) then
	538	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
	539	. mu0,fract,declin, zdtnirco2)
	540	endif
	541
	542	c Radiative flux from the sky absorbed by the surface (W.m-2)
	543	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	544	DO ig=1,ngrid
	545	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	546	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
	547	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
	548	ENDDO
	549
	550
	551	c Net atmospheric radiative heating rate (K.s-1)
	552	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	553	IF(callnlte) THEN
	554	CALL blendrad(ngrid, nlayer, pplay,
	555	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
	556	ELSE
	557	DO l=1,nlayer
	558	DO ig=1,ngrid
	559	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	560	& +zdtnirco2(ig,l)
	561	ENDDO
	562	ENDDO
	563	ENDIF
	564
	565	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	566
	567	c Transformation of the radiative tendencies:
	568	c -------------------------------------------
	569
	570	c Net radiative surface flux (W.m-2)
	571	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	572	c
	573	DO ig=1,ngrid
	574	zplanck(ig)=tsurf(ig)*tsurf(ig)
	575	zplanck(ig)=emis(ig)*
	576	$ stephanzplanck(ig)zplanck(ig)
	577	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
	578	ENDDO
	579
	580	DO l=1,nlayer
	581	DO ig=1,ngrid
	582	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	583	ENDDO
	584	ENDDO
	585
	586	ENDIF ! of IF (callrad)
	587
	588	c-----------------------------------------------------------------------
	589	c 3. Gravity wave and subgrid scale topography drag :
	590	c -------------------------------------------------
	591
	592
	593	IF(calllott)THEN
	594
	595	CALL calldrag_noro(ngrid,nlayer,ptimestep,
	596	& pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
	597
	598	DO l=1,nlayer
	599	DO ig=1,ngrid
	600	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
	601	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
	602	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
	603	ENDDO
	604	ENDDO
	605	ENDIF
	606
	607	c-----------------------------------------------------------------------
	608	c 4. Vertical diffusion (turbulent mixing):
	609	c -----------------------------------------
[161]	610
[38]	611	IF (calldifv) THEN
	612
	613	DO ig=1,ngrid
	614	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	615	ENDDO
	616
	617	CALL zerophys(ngrid*nlayer,zdum1)
	618	CALL zerophys(ngrid*nlayer,zdum2)
	619	do l=1,nlayer
	620	do ig=1,ngrid
	621	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	622	enddo
	623	enddo
[185]	624
[38]	625	c Calling vdif (Martian version WITH CO2 condensation)
	626	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
	627	$ ptimestep,capcal,lwrite,
	628	$ pplay,pplev,zzlay,zzlev,z0,
	629	$ pu,pv,zh,pq,tsurf,emis,qsurf,
	630	$ zdum1,zdum2,zdh,pdq,zflubid,
	631	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
	632	& zdqdif,zdqsdif)
	633
	634	DO l=1,nlayer
	635	DO ig=1,ngrid
	636	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	637	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	638	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	639
	640	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	641
	642	ENDDO
	643	ENDDO
	644
[185]	645	DO ig=1,ngrid
	646	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	647	ENDDO
[38]	648
	649	if (tracer) then
	650	DO iq=1, nq
	651	DO l=1,nlayer
	652	DO ig=1,ngrid
	653	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	654	ENDDO
	655	ENDDO
	656	ENDDO
	657	DO iq=1, nq
	658	DO ig=1,ngrid
	659	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	660	ENDDO
	661	ENDDO
	662	end if ! of if (tracer)
	663
	664	ELSE
	665	DO ig=1,ngrid
	666	zdtsurf(ig)=zdtsurf(ig)+
	667	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	668	ENDDO
	669	ENDIF ! of IF (calldifv)
	670
[161]	671	c-----------------------------------------------------------------------
	672	c TEST. Thermals :
	673	c HIGHLY EXPERIMENTAL, BEWARE !!
	674	c -----------------------------
	675
	676	if(calltherm) then
	677
[185]	678	call calltherm_interface(firstcall,
[161]	679	$ long,lati,zzlev,zzlay,
	680	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
[185]	681	$ pplay,pplev,pphi,zpopsk,
	682	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
	683	$ dtke_th,hfmax_th,wmax_th)
[161]	684
	685	DO l=1,nlayer
	686	DO ig=1,ngrid
	687	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
	688	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
	689	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
	690	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
	691	ENDDO
	692	ENDDO
	693
	694	DO ig=1,ngrid
	695	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
	696	ENDDO
	697
	698	if (tracer) then
	699	DO iq=1,nq
	700	DO l=1,nlayer
	701	DO ig=1,ngrid
	702	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
	703	ENDDO
	704	ENDDO
	705	ENDDO
	706	endif
[38]	707
[161]	708	else !of if calltherm
	709	lmax_th(:)=0
	710	end if
	711
[38]	712	c-----------------------------------------------------------------------
	713	c 5. Dry convective adjustment:
	714	c -----------------------------
	715
	716	IF(calladj) THEN
	717
	718	DO l=1,nlayer
	719	DO ig=1,ngrid
	720	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	721	ENDDO
	722	ENDDO
	723	CALL zerophys(ngrid*nlayer,zduadj)
	724	CALL zerophys(ngrid*nlayer,zdvadj)
	725	CALL zerophys(ngrid*nlayer,zdhadj)
	726
	727	CALL convadj(ngrid,nlayer,nq,ptimestep,
[161]	728	$ pplay,pplev,zpopsk,lmax_th,
[38]	729	$ pu,pv,zh,pq,
	730	$ pdu,pdv,zdh,pdq,
	731	$ zduadj,zdvadj,zdhadj,
	732	$ zdqadj)
	733
[161]	734
[38]	735	DO l=1,nlayer
	736	DO ig=1,ngrid
	737	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	738	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	739	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	740
	741	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	742	ENDDO
	743	ENDDO
	744
	745	if(tracer) then
	746	DO iq=1, nq
	747	DO l=1,nlayer
	748	DO ig=1,ngrid
	749	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	750	ENDDO
	751	ENDDO
	752	ENDDO
	753	end if
	754	ENDIF ! of IF(calladj)
	755
	756	c-----------------------------------------------------------------------
	757	c 6. Carbon dioxide condensation-sublimation:
	758	c -------------------------------------------
	759
	760	IF (callcond) THEN
	761	CALL newcondens(ngrid,nlayer,nq,ptimestep,
	762	$ capcal,pplay,pplev,tsurf,pt,
	763	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	764	$ co2ice,albedo,emis,
	765	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
[185]	766	$ fluxsurf_sw,zls)
[38]	767
	768	DO l=1,nlayer
	769	DO ig=1,ngrid
	770	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	771	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	772	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	773	ENDDO
	774	ENDDO
	775	DO ig=1,ngrid
	776	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	777	ENDDO
	778
	779	IF (tracer) THEN
	780	DO iq=1, nq
	781	DO l=1,nlayer
	782	DO ig=1,ngrid
	783	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	784	ENDDO
	785	ENDDO
	786	ENDDO
	787	ENDIF ! of IF (tracer)
	788
	789	ENDIF ! of IF (callcond)
	790
	791	c-----------------------------------------------------------------------
	792	c 7. Specific parameterizations for tracers
	793	c: -----------------------------------------
	794
	795	if (tracer) then
	796
	797	c 7a. Water and ice
	798	c ---------------
	799
	800	c ---------------------------------------
	801	c Water ice condensation in the atmosphere
	802	c ----------------------------------------
	803	IF (water) THEN
	804
	805	call watercloud(ngrid,nlayer,ptimestep,
	806	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
	807	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
	808	& nq,naerkind,tau,
	809	& ccn,rdust,rice,nuice)
	810	if (activice) then
	811	c Temperature variation due to latent heat release
	812	DO l=1,nlayer
	813	DO ig=1,ngrid
	814	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
	815	ENDDO
	816	ENDDO
	817	endif
	818
	819	! increment water vapour and ice atmospheric tracers tendencies
	820	IF (water) THEN
	821	DO l=1,nlayer
	822	DO ig=1,ngrid
	823	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	824	& zdqcloud(ig,l,igcm_h2o_vap)
	825	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	826	& zdqcloud(ig,l,igcm_h2o_ice)
	827	ENDDO
	828	ENDDO
	829	ENDIF ! of IF (water) THEN
	830	! Increment water ice surface tracer tendency
	831	DO ig=1,ngrid
	832	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
	833	& zdqscloud(ig,igcm_h2o_ice)
	834	ENDDO
	835
	836	END IF ! of IF (water)
	837
	838	c 7b. Chemical species
	839	c ------------------
	840
	841	c --------------
	842	c photochemistry :
	843	c --------------
	844	IF (photochem .or. thermochem) then
	845	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
	846	& zzlay,zday,pq,pdq,rice,
	847	& zdqchim,zdqschim,zdqcloud,zdqscloud)
	848	!NB: Photochemistry includes condensation of H2O2
	849
	850	! increment values of tracers:
	851	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	852	! tracers is zero anyways
	853	DO l=1,nlayer
	854	DO ig=1,ngrid
	855	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	856	ENDDO
	857	ENDDO
	858	ENDDO ! of DO iq=1,nq
	859	! add condensation tendency for H2O2
	860	if (igcm_h2o2.ne.0) then
	861	DO l=1,nlayer
	862	DO ig=1,ngrid
	863	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	864	& +zdqcloud(ig,l,igcm_h2o2)
	865	ENDDO
	866	ENDDO
	867	endif
	868
	869	! increment surface values of tracers:
	870	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	871	! tracers is zero anyways
	872	DO ig=1,ngrid
	873	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	874	ENDDO
	875	ENDDO ! of DO iq=1,nq
	876	! add condensation tendency for H2O2
	877	if (igcm_h2o2.ne.0) then
	878	DO ig=1,ngrid
	879	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	880	& +zdqscloud(ig,igcm_h2o2)
	881	ENDDO
	882	endif
	883
	884	END IF ! of IF (photochem.or.thermochem)
	885
	886	c 7c. Aerosol particles
	887	c -------------------
	888
	889	c ----------
	890	c Dust devil :
	891	c ----------
	892	IF(callddevil) then
	893	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
	894	& zdqdev,zdqsdev)
	895
	896	if (dustbin.ge.1) then
	897	do iq=1,nq
	898	DO l=1,nlayer
	899	DO ig=1,ngrid
	900	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	901	ENDDO
	902	ENDDO
	903	enddo
	904	do iq=1,nq
	905	DO ig=1,ngrid
	906	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
	907	ENDDO
	908	enddo
	909	endif ! of if (dustbin.ge.1)
	910
	911	END IF ! of IF (callddevil)
	912
	913	c -------------
	914	c Sedimentation : acts also on water ice
	915	c -------------
	916	IF (sedimentation) THEN
	917	!call zerophys(ngridnlayernq, zdqsed)
	918	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	919	!call zerophys(ngrid*nq, zdqssed)
	920	zdqssed(1:ngrid,1:nq)=0
	921
	922	call callsedim(ngrid,nlayer, ptimestep,
	923	& pplev,zzlev, pt, rdust, rice,
	924	& pq, pdq, zdqsed, zdqssed,nq)
	925
	926	DO iq=1, nq
	927	DO l=1,nlayer
	928	DO ig=1,ngrid
	929	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	930	ENDDO
	931	ENDDO
	932	ENDDO
	933	DO iq=1, nq
	934	DO ig=1,ngrid
	935	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	936	ENDDO
	937	ENDDO
	938	END IF ! of IF (sedimentation)
	939
	940	c 7d. Updates
	941	c ---------
	942
	943	DO iq=1, nq
	944	DO ig=1,ngrid
	945
	946	c ---------------------------------
	947	c Updating tracer budget on surface
	948	c ---------------------------------
	949	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	950
	951	ENDDO ! (ig)
	952	ENDDO ! (iq)
	953
	954	endif ! of if (tracer)
	955
	956
	957	c-----------------------------------------------------------------------
	958	c 8. THERMOSPHERE CALCULATION
	959	c-----------------------------------------------------------------------
	960
	961	if (callthermos) then
	962	call thermosphere(pplev,pplay,dist_sol,
	963	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	964	& pt,pq,pu,pv,pdt,pdq,
	965	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	966
	967	DO l=1,nlayer
	968	DO ig=1,ngrid
	969	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	970	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	971	& +zdteuv(ig,l)
	972	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	973	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	974	DO iq=1, nq
	975	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	976	ENDDO
	977	ENDDO
	978	ENDDO
	979
	980	endif ! of if (callthermos)
	981
	982	c-----------------------------------------------------------------------
	983	c 9. Surface and sub-surface soil temperature
	984	c-----------------------------------------------------------------------
	985	c
	986	c
	987	c 9.1 Increment Surface temperature:
	988	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	989
	990	DO ig=1,ngrid
	991	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	992	ENDDO
	993
	994	c Prescribe a cold trap at south pole (except at high obliquity !!)
	995	c Temperature at the surface is set there to be the temperature
	996	c corresponding to equilibrium temperature between phases of CO2
	997
	998	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
	999	if (caps.and.(obliquit.lt.27.)) then
	1000	! NB: Updated surface pressure, at grid point 'ngrid', is
	1001	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
	1002	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	1003	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	1004	endif
	1005	c -------------------------------------------------------------
	1006	c Change of surface albedo (set to 0.4) in case of ground frost
	1007	c everywhere except on the north permanent cap and in regions
	1008	c covered by dry ice.
	1009	c ALWAYS PLACE these lines after newcondens !!!
	1010	c -------------------------------------------------------------
	1011	do ig=1,ngrid
	1012	if ((co2ice(ig).eq.0).and.
	1013	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
	1014	albedo(ig,1) = alb_surfice
	1015	albedo(ig,2) = alb_surfice
	1016	endif
	1017	enddo ! of do ig=1,ngrid
	1018	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
	1019
	1020	c
	1021	c 9.2 Compute soil temperatures and subsurface heat flux:
	1022	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1023	IF (callsoil) THEN
	1024	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
	1025	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1026	ENDIF
	1027
	1028	c-----------------------------------------------------------------------
	1029	c 10. Write output files
	1030	c ----------------------
	1031
	1032	c -------------------------------
	1033	c Dynamical fields incrementation
	1034	c -------------------------------
	1035	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1036	! temperature, zonal and meridional wind
	1037	DO l=1,nlayer
	1038	DO ig=1,ngrid
	1039	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1040	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1041	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1042	ENDDO
	1043	ENDDO
	1044
	1045	! tracers
	1046	DO iq=1, nq
	1047	DO l=1,nlayer
	1048	DO ig=1,ngrid
	1049	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1050	ENDDO
	1051	ENDDO
	1052	ENDDO
	1053
	1054	! surface pressure
	1055	DO ig=1,ngrid
	1056	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1057	ENDDO
	1058
	1059	! pressure
	1060	DO l=1,nlayer
	1061	DO ig=1,ngrid
	1062	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1063	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1064	ENDDO
	1065	ENDDO
	1066
	1067	! Density
	1068	DO l=1,nlayer
	1069	DO ig=1,ngrid
	1070	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1071	ENDDO
	1072	ENDDO
	1073
	1074	c Compute surface stress : (NB: z0 is a common in planete.h)
	1075	c DO ig=1,ngrid
	1076	c cd = (0.4/log(zzlay(ig,1)/z0))**2
	1077	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1078	c ENDDO
	1079
	1080	c Sum of fluxes in solar spectral bands (for output only)
	1081	DO ig=1,ngrid
	1082	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1083	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
	1084	ENDDO
	1085	c ***** TEST ****************************************************
	1086	ztim1 = 999
	1087	DO l=1,nlayer
	1088	DO ig=1,ngrid
	1089	if (pt(ig,l).lt.ztim1) then
	1090	ztim1 = pt(ig,l)
	1091	igmin = ig
	1092	lmin = l
	1093	end if
	1094	ENDDO
	1095	ENDDO
	1096	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
	1097	write(,) 'PHYSIQ: stability WARNING :'
	1098	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1099	& 'ig l =', igmin, lmin
	1100	end if
	1101	c *******************************************************************
	1102
	1103	c ---------------------
	1104	c Outputs to the screen
	1105	c ---------------------
	1106
	1107	IF (lwrite) THEN
	1108	PRINT*,'Global diagnostics for the physics'
	1109	PRINT,'Variables and their increments x and dx/dt dt'
	1110	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1111	WRITE(*,'(2f10.5,2f15.5)')
	1112	s tsurf(igout),zdtsurf(igout)*ptimestep,
	1113	s pplev(igout,1),pdpsrf(igout)*ptimestep
	1114	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1115	WRITE(*,'(i4,6f10.5)') (l,
	1116	s pu(igout,l),pdu(igout,l)*ptimestep,
	1117	s pv(igout,l),pdv(igout,l)*ptimestep,
	1118	s pt(igout,l),pdt(igout,l)*ptimestep,
	1119	s l=1,nlayer)
	1120	ENDIF ! of IF (lwrite)
	1121
	1122	IF (ngrid.NE.1) THEN
	1123	print,'Ls =',zls180./pi,
	1124	& ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
	1125	& ' tau(Viking1) =',tau(ig_vl1,1)
	1126
	1127
	1128	c -------------------------------------------------------------------
	1129	c Writing NetCDF file "RESTARTFI" at the end of the run
	1130	c -------------------------------------------------------------------
	1131	c Note: 'restartfi' is stored just before dynamics are stored
	1132	c in 'restart'. Between now and the writting of 'restart',
	1133	c there will have been the itau=itau+1 instruction and
	1134	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1135	c thus we store for time=time+dtvr
	1136
	1137	IF(lastcall) THEN
	1138	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1139	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1140	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1141	. ptimestep,pday,
	1142	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
	1143	. area,albedodat,inertiedat,zmea,zstd,zsig,
	1144	. zgam,zthe)
	1145	ENDIF
	1146
	1147	c -------------------------------------------------------------------
	1148	c Calculation of diagnostic variables written in both stats and
	1149	c diagfi files
	1150	c -------------------------------------------------------------------
	1151
	1152	if (tracer) then
	1153	if (water) then
	1154
	1155	call zerophys(ngrid,mtot)
	1156	call zerophys(ngrid,icetot)
	1157	call zerophys(ngrid,rave)
	1158	call zerophys(ngrid,tauTES)
	1159	do ig=1,ngrid
	1160	do l=1,nlayermx
	1161	mtot(ig) = mtot(ig) +
	1162	& zq(ig,l,igcm_h2o_vap) *
	1163	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1164	icetot(ig) = icetot(ig) +
	1165	& zq(ig,l,igcm_h2o_ice) *
	1166	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1167	rave(ig) = rave(ig) +
	1168	& zq(ig,l,igcm_h2o_ice) *
	1169	& (pplev(ig,l) - pplev(ig,l+1)) / g *
	1170	& rice(ig,l) * (1.+nuice_ref)
	1171	c Computing abs optical depth at 825 cm-1 in each
	1172	c layer to simulate NEW TES retrieval
	1173	Qabsice = min(
	1174	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1175	& )
	1176	opTES(ig,l)= 0.75 * Qabsice *
	1177	& zq(ig,l,igcm_h2o_ice) *
	1178	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1179	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1180	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1181	enddo
	1182	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
	1183	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1184	enddo
	1185
	1186	endif ! of if (water)
	1187	endif ! of if (tracer)
	1188
	1189	c -----------------------------------------------------------------
	1190	c WSTATS: Saving statistics
	1191	c -----------------------------------------------------------------
	1192	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1193	c which can later be used to make the statistic files of the run:
	1194	c "stats") only possible in 3D runs !
	1195
	1196	IF (callstats) THEN
	1197
	1198	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1199	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1200	call wstats(ngrid,"co2ice","CO2 ice cover",
	1201	& "kg.m-2",2,co2ice)
	1202	call wstats(ngrid,"fluxsurf_lw",
	1203	& "Thermal IR radiative flux to surface","W.m-2",2,
	1204	& fluxsurf_lw)
	1205	call wstats(ngrid,"fluxsurf_sw",
	1206	& "Solar radiative flux to surface","W.m-2",2,
	1207	& fluxsurf_sw_tot)
	1208	call wstats(ngrid,"fluxtop_lw",
	1209	& "Thermal IR radiative flux to space","W.m-2",2,
	1210	& fluxtop_lw)
	1211	call wstats(ngrid,"fluxtop_sw",
	1212	& "Solar radiative flux to space","W.m-2",2,
	1213	& fluxtop_sw_tot)
	1214	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1215	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1216	call wstats(ngrid,"v","Meridional (North-South) wind",
	1217	& "m.s-1",3,zv)
	1218	call wstats(ngrid,"w","Vertical (down-up) wind",
	1219	& "m.s-1",3,pw)
	1220	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
	1221	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
	1222	c call wstats(ngrid,"q2",
	1223	c & "Boundary layer eddy kinetic energy",
	1224	c & "m2.s-2",3,q2)
	1225	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1226	c & emis)
	1227	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1228	c & 2,zstress)
	1229	c call wstats(ngrid,"sw_htrt","sw heat.rate",
	1230	c & "W.m-2",3,zdtsw)
	1231	c call wstats(ngrid,"lw_htrt","lw heat.rate",
	1232	c & "W.m-2",3,zdtlw)
	1233
	1234	if (tracer) then
	1235	if (water) then
	1236	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1237	& *mugaz/mmol(igcm_h2o_vap)
	1238	call wstats(ngrid,"vmr_h2ovapor",
	1239	& "H2O vapor volume mixing ratio","mol/mol",
	1240	& 3,vmr)
	1241	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1242	& *mugaz/mmol(igcm_h2o_ice)
	1243	call wstats(ngrid,"vmr_h2oice",
	1244	& "H2O ice volume mixing ratio","mol/mol",
	1245	& 3,vmr)
	1246	call wstats(ngrid,"h2o_ice_s",
	1247	& "surface h2o_ice","kg/m2",
	1248	& 2,qsurf(1,igcm_h2o_ice))
	1249
	1250	call wstats(ngrid,"mtot",
	1251	& "total mass of water vapor","kg/m2",
	1252	& 2,mtot)
	1253	call wstats(ngrid,"icetot",
	1254	& "total mass of water ice","kg/m2",
	1255	& 2,icetot)
	1256	call wstats(ngrid,"reffice",
	1257	& "Mean reff","m",
	1258	& 2,rave)
	1259	c call wstats(ngrid,"rice",
	1260	c & "Ice particle size","m",
	1261	c & 3,rice)
	1262	c If activice is true, tauTES is computed in aeropacity.F;
	1263	if (.not.activice) then
	1264	call wstats(ngrid,"tauTESap",
	1265	& "tau abs 825 cm-1","",
	1266	& 2,tauTES)
	1267	endif
	1268
	1269	endif ! of if (water)
	1270
	1271	if (thermochem.or.photochem) then
	1272	do iq=1,nq
	1273	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
	1274	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
	1275	. (noms(iq).eq."h2").or.
	1276	. (noms(iq).eq."o3")) then
	1277	do l=1,nlayer
	1278	do ig=1,ngrid
	1279	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
	1280	end do
	1281	end do
	1282	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1283	. "Volume mixing ratio","mol/mol",3,vmr)
	1284	endif
	1285	enddo
	1286	endif ! of if (thermochem.or.photochem)
	1287
	1288	endif ! of if (tracer)
	1289
	1290	IF(lastcall) THEN
	1291	write (,) "Writing stats..."
	1292	call mkstats(ierr)
	1293	ENDIF
	1294
	1295	ENDIF !if callstats
	1296
	1297	c (Store EOF for Mars Climate database software)
	1298	IF (calleofdump) THEN
	1299	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1300	ENDIF
	1301
	1302	c ==========================================================
	1303	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1304	c any variable for diagnostic (output with period
	1305	c "ecritphy", set in "run.def")
	1306	c ==========================================================
	1307	c WRITEDIAGFI can ALSO be called from any other subroutines
	1308	c for any variables !!
	1309	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1310	c & emis)
[164]	1311	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
	1312	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
[38]	1313	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1314	& tsurf)
	1315	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
[164]	1316	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1317	& co2ice)
[38]	1318	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1319	c & "K",2,zt(1,7))
	1320	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1321	& fluxsurf_lw)
	1322	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1323	& fluxsurf_sw_tot)
	1324	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1325	& fluxtop_lw)
	1326	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1327	& fluxtop_sw_tot)
	1328	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
[161]	1329	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1330	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1331	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
[38]	1332	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
	1333	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
[164]	1334	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
[38]	1335	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
	1336	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1337	c & zstress)
	1338	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	1339	c & 'w.m-2',3,zdtsw)
	1340	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	1341	c & 'w.m-2',3,zdtlw)
[161]	1342	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1343	c & 'tau abs 825 cm-1',
	1344	c & '',2,tauTES)
[38]	1345
[161]	1346
[38]	1347	c ----------------------------------------------------------
	1348	c Outputs of the CO2 cycle
	1349	c ----------------------------------------------------------
	1350
	1351	if (tracer.and.(igcm_co2.ne.0)) then
	1352	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	1353	! & "kg/kg",2,zq(1,1,igcm_co2))
	1354	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	1355	! & "kg/kg",3,zq(1,1,igcm_co2))
	1356
	1357	! Compute co2 column
	1358	call zerophys(ngrid,co2col)
	1359	do l=1,nlayermx
	1360	do ig=1,ngrid
	1361	co2col(ig)=co2col(ig)+
	1362	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
	1363	enddo
	1364	enddo
	1365	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	1366	& co2col)
	1367	endif ! of if (tracer.and.(igcm_co2.ne.0))
	1368
	1369	c ----------------------------------------------------------
	1370	c Outputs of the water cycle
	1371	c ----------------------------------------------------------
	1372	if (tracer) then
	1373	if (water) then
	1374
	1375	CALL WRITEDIAGFI(ngridmx,'mtot',
	1376	& 'total mass of water vapor',
	1377	& 'kg/m2',2,mtot)
	1378	CALL WRITEDIAGFI(ngridmx,'icetot',
	1379	& 'total mass of water ice',
	1380	& 'kg/m2',2,icetot)
	1381	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1382	c & *mugaz/mmol(igcm_h2o_ice)
	1383	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
	1384	c & 'mol/mol',3,vmr)
	1385	CALL WRITEDIAGFI(ngridmx,'reffice',
	1386	& 'Mean reff',
	1387	& 'm',2,rave)
	1388	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
	1389	c & 'm',3,rice)
	1390	c If activice is true, tauTES is computed in aeropacity.F;
	1391	if (.not.activice) then
	1392	CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1393	& 'tau abs 825 cm-1',
	1394	& '',2,tauTES)
	1395	endif
	1396	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
	1397	& 'surface h2o_ice',
	1398	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
	1399	endif !(water)
	1400
	1401
	1402	if (water.and..not.photochem) then
	1403	iq=nq
	1404	c write(str2(1:2),'(i2.2)') iq
	1405	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
	1406	c & 'kg.m-2',2,zdqscloud(1,iq))
	1407	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
	1408	c & 'kg/kg',3,zdqchim(1,1,iq))
	1409	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
	1410	c & 'kg/kg',3,zdqdif(1,1,iq))
	1411	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
	1412	c & 'kg/kg',3,zdqadj(1,1,iq))
	1413	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
	1414	c & 'kg/kg',3,zdqc(1,1,iq))
	1415	endif !(water.and..not.photochem)
	1416	endif
	1417
	1418	c ----------------------------------------------------------
	1419	c Outputs of the dust cycle
	1420	c ----------------------------------------------------------
	1421
	1422	c call WRITEDIAGFI(ngridmx,'tauref',
	1423	c & 'Dust ref opt depth','NU',2,tauref)
	1424
	1425	if (tracer.and.(dustbin.ne.0)) then
	1426	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
	1427	if (doubleq) then
	1428	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
	1429	c & 'kg.m-2',2,qsurf(1,1))
	1430	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
	1431	c & 'N.m-2',2,qsurf(1,2))
	1432	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
	1433	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
	1434	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
	1435	c & 'kg.m-2.s-1',2,zdqssed(1,1))
	1436	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
	1437	& 'm',3,rdust*ref_r0)
	1438	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	1439	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
	1440	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
	1441	c & 'part/kg',3,pq(1,1,igcm_dust_number))
	1442	else
	1443	do iq=1,dustbin
	1444	write(str2(1:2),'(i2.2)') iq
	1445	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
	1446	& 'kg/kg',3,zq(1,1,iq))
	1447	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
	1448	& 'kg.m-2',2,qsurf(1,iq))
	1449	end do
	1450	endif ! (doubleq)
	1451	c if (submicron) then
	1452	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
	1453	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
	1454	c endif ! (submicron)
	1455	end if ! (tracer.and.(dustbin.ne.0))
	1456
	1457	c ----------------------------------------------------------
[161]	1458	c Outputs of thermals
	1459	c ----------------------------------------------------------
[173]	1460	if (calltherm) then
[161]	1461
	1462	! call WRITEDIAGFI(ngrid,'dtke',
	1463	! & 'tendance tke thermiques','m2/s2',
	1464	! & 3,dtke_th)
	1465	! call WRITEDIAGFI(ngrid,'d_u_ajs',
	1466	! & 'tendance u thermiques','m/s',
	1467	! & 3,pdu_th*ptimestep)
	1468	! call WRITEDIAGFI(ngrid,'d_v_ajs',
	1469	! & 'tendance v thermiques','m/s',
	1470	! & 3,pdv_th*ptimestep)
	1471	! if (tracer) then
	1472	! if (nq .eq. 2) then
	1473	! call WRITEDIAGFI(ngrid,'deltaq_th',
	1474	! & 'delta q thermiques','kg/kg',
	1475	! & 3,ptimestep*pdq_th(:,:,2))
	1476	! endif
	1477	! endif
	1478
[173]	1479	lmax_th_out(:)=real(lmax_th(:))
[161]	1480
[173]	1481	call WRITEDIAGFI(ngridmx,'lmax_th',
	1482	& 'hauteur du thermique','K',
	1483	& 2,lmax_th_out)
	1484	call WRITEDIAGFI(ngridmx,'hfmax_th',
	1485	& 'maximum TH heat flux','K.m/s',
	1486	& 2,hfmax_th)
	1487	call WRITEDIAGFI(ngridmx,'wmax_th',
	1488	& 'maximum TH vertical velocity','m/s',
	1489	& 2,wmax_th)
	1490
	1491	endif
	1492
[161]	1493	c ----------------------------------------------------------
[38]	1494	c Output in netcdf file "diagsoil.nc" for subterranean
	1495	c variables (output every "ecritphy", as for writediagfi)
	1496	c ----------------------------------------------------------
	1497
	1498	! Write soil temperature
	1499	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	1500	! & 3,tsoil)
	1501	! Write surface temperature
	1502	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	1503	! & 2,tsurf)
	1504
	1505	c ==========================================================
	1506	c END OF WRITEDIAGFI
	1507	c ==========================================================
	1508
	1509	ELSE ! if(ngrid.eq.1)
	1510
	1511	print,'Ls =',zls180./pi,
	1512	& ' tauref(700 Pa) =',tauref
	1513	c ----------------------------------------------------------------------
	1514	c Output in grads file "g1d" (ONLY when using testphys1d)
	1515	c (output at every X physical timestep)
	1516	c ----------------------------------------------------------------------
	1517	c
	1518	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	1519	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	1520	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	1521
	1522	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	1523	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	1524	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	1525	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	1526
[161]	1527	! THERMALS STUFF 1D
[173]	1528
[161]	1529	if(calltherm) then
	1530
	1531	lmax_th_out(:)=real(lmax_th(:))
	1532
	1533	call WRITEDIAGFI(ngridmx,'lmax_th',
[173]	1534	& 'hauteur du thermique','point',
[161]	1535	& 0,lmax_th_out)
[173]	1536	call WRITEDIAGFI(ngridmx,'hfmax_th',
	1537	& 'maximum TH heat flux','K.m/s',
	1538	& 0,hfmax_th)
	1539	call WRITEDIAGFI(ngridmx,'wmax_th',
	1540	& 'maximum TH vertical velocity','m/s',
	1541	& 0,wmax_th)
[161]	1542
	1543
	1544	co2col(:)=0.
	1545	dummycol(:)=0.
	1546	if (tracer) then
	1547	do l=1,nlayermx
	1548	do ig=1,ngrid
	1549	co2col(ig)=co2col(ig)+
	1550	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
	1551	if (nqmx .gt. 1) then
[164]	1552	iq=2 ! to avoid out-of-bounds spotting by picky compilers
	1553	! when gcm is compiled with only one tracer
[161]	1554	dummycol(ig)=dummycol(ig)+
[164]	1555	& zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g
[161]	1556	endif
	1557	enddo
	1558	enddo
	1559
	1560	end if
	1561	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
	1562	& ,'kg/m-2',0,co2col)
	1563	call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass' &
	1564	& ,'kg/m-2',0,dummycol)
	1565	endif
	1566	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
	1567	& ,'m/s',1,pw)
	1568	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
	1569	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
	1570	& tsurf)
	1571
	1572	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
	1573	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
[38]	1574	! or output in diagfi.nc (for testphys1d)
	1575	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
	1576	call WRITEDIAGFI(ngridmx,'temp','Temperature',
	1577	& 'K',1,zt)
	1578
	1579	if(tracer) then
	1580	c CALL writeg1d(ngrid,1,tau,'tau','SI')
	1581	do iq=1,nq
	1582	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	1583	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
	1584	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	1585	end do
	1586	end if
	1587
	1588	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
	1589
	1590	do l=2,nlayer-1
	1591	tmean=zt(1,l)
	1592	if(zt(1,l).ne.zt(1,l-1))
	1593	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	1594	zlocal(l)= zlocal(l-1)
	1595	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
	1596	enddo
	1597	zlocal(nlayer)= zlocal(nlayer-1)-
	1598	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
	1599	& rnew(1,nlayer)*tmean/g
	1600
	1601	END IF ! if(ngrid.ne.1)
	1602
	1603	icount=icount+1
	1604	RETURN
	1605	END

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

Download in other formats:

Original Format