Context Navigation

physiq.F @ 288

Last change on this file since 288 was 288, checked in by aslmd, 13 years ago

LMDZ.MARS

09/09/11 == AS

--> Fixed a problem with .eq. used with booleans in physiq.F
[some good picky compilers complain about this]
--> Added a warning in inifis about using
callrichsl = .false.
calltherm = .true.
which is not recommended. The new surface layer has been built to go
with the new mixing layer scheme (thermals). And anyway it is a much
better scheme than the previous one.

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

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

Download in other formats:

Original Format