Context Navigation

physiq.F @ 258

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

--- AC 03/08/2011 ---

M 255 libf/phymars/physiq.F
<> Modified to interface new SL parametrization

M 255 libf/phymars/vdif_cd.F
<> New SL parametrization based on a bulk Richardson Monin-Obukhov theory formulation

Stability functions are taken from D.E. England et al. (95)
Similarity functions coefficients based on Dyer and Hicks (70).
Includes thermal roughness length computation, heat and momentum drag coefficient computation
Can be used to output hydrodynamic-related SL quantities (bulk Richardson, turbulent Prandtl number estimation, Reynolds number)

M 255 libf/phymars/thermcell_dqupdown.F90
<> Minor modification to suit picky compilers

M 255 libf/phymars/vdifc.F
<> Now takes into account sub-grid gustiness, evaluated from thermals activity (it's proxy being the maximum vertical velocity)

M 255 libf/phymars/meso_inc/meso_inc_les.F
<> Minor modification : u* is now taken from the new vdifc and not recomputed from a simple law

File size: 60.2 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 "watercap.h"
	131	#include "param.h"
	132	#include "param_v3.h"
	133	#include "conc.h"
	134
	135	#include "netcdf.inc"
	136
[234]	137	#include "slope.h"
	138
[86]	139	#ifdef MESOSCALE
[42]	140	#include "wrf_output_2d.h"
	141	#include "wrf_output_3d.h"
[226]	142	#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
	143	#include "meso_inc/meso_inc_var.F"
[86]	144	#endif
[42]	145
	146	c Arguments :
	147	c -----------
	148
	149	c inputs:
	150	c -------
	151	INTEGER ngrid,nlayer,nq
	152	REAL ptimestep
	153	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	154	REAL pphi(ngridmx,nlayer)
	155	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	156	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	157	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
	158	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	159	LOGICAL firstcall,lastcall
[226]	160
[42]	161	REAL pday
	162	REAL ptime
	163	logical tracerdyn
	164
	165	c outputs:
	166	c --------
	167	c physical tendencies
	168	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	169	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	170	REAL pdpsrf(ngridmx) ! surface pressure tendency
	171
	172
	173	c Local saved variables:
	174	c ----------------------
	175	c aerosol (dust or ice) extinction optical depth at reference wavelength
	176	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
	177	c aerosol optical properties :
	178	REAL aerosol(ngridmx,nlayermx,naerkind)
	179
	180	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	181	INTEGER icount ! counter of calls to physiq during the run.
	182	REAL tsurf(ngridmx) ! Surface temperature (K)
	183	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	184	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
	185	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
	186	REAL emis(ngridmx) ! Thermal IR surface emissivity
	187	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	188	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	189	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	190	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	191	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
	192	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	193	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	194
	195	c Variables used by the water ice microphysical scheme:
	196	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
	197	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
	198	! of the size distribution
	199	c Albedo of deposited surface ice
[82]	200	!!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
[162]	201	REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
[42]	202
[234]	203	c Variables used by the slope model
	204	REAL sl_ls, sl_lct, sl_lat
	205	REAL sl_tau, sl_alb, sl_the, sl_psi
	206	REAL sl_fl0, sl_flu
	207	REAL sl_ra, sl_di0
	208	REAL sky
	209
[42]	210	SAVE day_ini, icount
	211	SAVE aerosol, tsurf,tsoil
	212	SAVE co2ice,albedo,emis, q2
	213	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
	214
	215	REAL stephan
	216	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	217	SAVE stephan
	218
	219	c Local variables :
	220	c -----------------
	221
	222	REAL CBRT
	223	EXTERNAL CBRT
	224
	225	CHARACTER*80 fichier
	226	INTEGER l,ig,ierr,igout,iq,i, tapphys
	227
	228	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
	229	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
	230	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
	231	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
	232	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
	233	! (used if active=F)
	234	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	235	REAL zls ! solar longitude (rad)
	236	REAL zday ! date (time since Ls=0, in martian days)
	237	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	238	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
	239	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
	240
	241	c Tendancies due to various processes:
	242	REAL dqsurf(ngridmx,nqmx)
	243	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	244	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
	245	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
	246	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
	247	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
	248	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
	249	REAL zdtsurf(ngridmx) ! (K/s)
	250	REAL zdtcloud(ngridmx,nlayermx)
	251	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	252	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	253	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	254	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	255	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	256	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	257	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	258	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	259
	260	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
	261	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
	262	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
	263	REAL zdqadj(ngridmx,nlayermx,nqmx)
	264	REAL zdqc(ngridmx,nlayermx,nqmx)
	265	REAL zdqcloud(ngridmx,nlayermx,nqmx)
	266	REAL zdqscloud(ngridmx,nqmx)
	267	REAL zdqchim(ngridmx,nlayermx,nqmx)
	268	REAL zdqschim(ngridmx,nqmx)
	269
	270	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	271	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	272	REAL zdumolvis(ngridmx,nlayermx)
	273	REAL zdvmolvis(ngridmx,nlayermx)
	274	real zdqmoldiff(ngridmx,nlayermx,nqmx)
	275
	276	c Local variable for local intermediate calcul:
	277	REAL zflubid(ngridmx)
	278	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	279	REAL zdum1(ngridmx,nlayermx)
	280	REAL zdum2(ngridmx,nlayermx)
	281	REAL ztim1,ztim2,ztim3, z1,z2
	282	REAL ztime_fin
	283	REAL zdh(ngridmx,nlayermx)
	284	INTEGER length
	285	PARAMETER (length=100)
	286
	287	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	288	c -----------------------------------------------------------------------------
	289	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	290	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	291	REAL zq(ngridmx,nlayermx,nqmx)
	292	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	293	character*2 str2
	294	character*5 str5
	295	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
	296	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
	297	! (particules kg-1)
[162]	298	SAVE ccn !! in case iradia != 1
[42]	299	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
	300	real qtot1,qtot2 ! total aerosol mass
	301	integer igmin, lmin
	302	logical tdiag
	303
	304	real co2col(ngridmx) ! CO2 column
	305	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
	306	REAL zstress(ngrid), cd
	307	real hco2(nqmx),tmean, zlocal(nlayermx)
	308	real rho(ngridmx,nlayermx) ! density
	309	real vmr(ngridmx,nlayermx) ! volume mixing ratio
	310	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
	311	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
	312	REAL rave(ngridmx) ! Mean water ice effective radius (m)
	313	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
	314	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
	315	REAL Qabsice ! Water ice absorption coefficient
	316
	317
	318	REAL time_phys
	319
[162]	320	c Variables from thermal
	321
[185]	322	REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
[256]	323	REAL, SAVE :: wmax_th(ngridmx)
[234]	324	REAL, SAVE :: hfmax_th(ngridmx)
[185]	325	REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
	326	REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
	327	INTEGER lmax_th(ngridmx)
	328	REAL dtke_th(ngridmx,nlayermx+1)
[256]	329	REAL zcdv(ngridmx)
[42]	330	c=======================================================================
	331
	332	c 1. Initialisation:
	333	c -----------------
	334
	335	c 1.1 Initialisation only at first call
	336	c ---------------------------------------
	337	IF (firstcall) THEN
	338
	339	c variables set to 0
	340	c ~~~~~~~~~~~~~~~~~~
	341	call zerophys(ngridnlayernaerkind,aerosol)
	342	call zerophys(ngrid*nlayer,dtrad)
	343	call zerophys(ngrid,fluxrad)
	344
	345	c read startfi
	346	c ~~~~~~~~~~~~
[226]	347	#ifndef MESOSCALE
	348	! Read netcdf initial physical parameters.
	349	CALL phyetat0 ("startfi.nc",0,0,
	350	& nsoilmx,nq,
	351	& day_ini,time_phys,
	352	& tsurf,tsoil,emis,q2,qsurf,co2ice)
	353	#else
	354	#include "meso_inc/meso_inc_ini.F"
	355	#endif
[185]	356
[42]	357	if (pday.ne.day_ini) then
	358	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	359	& "physics and dynamics"
	360	write(,) "dynamics day: ",pday
	361	write(,) "physics day: ",day_ini
	362	stop
	363	endif
	364
	365	write (,) 'In physiq day_ini =', day_ini
	366
	367	c Initialize albedo and orbital calculation
	368	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	369	CALL surfini(ngrid,co2ice,qsurf,albedo)
	370	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	371
	372	c initialize soil
	373	c ~~~~~~~~~~~~~~~
	374	IF (callsoil) THEN
	375	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	376	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	377	ELSE
	378	PRINT*,
	379	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	380	DO ig=1,ngrid
	381	capcal(ig)=1.e5
	382	fluxgrd(ig)=0.
	383	ENDDO
	384	ENDIF
	385	icount=1
	386
	387
	388	c initialize tracers
	389	c ~~~~~~~~~~~~~~~~~~
	390	tracerdyn=tracer
	391	IF (tracer) THEN
	392	CALL initracer(qsurf,co2ice)
	393	ENDIF ! end tracer
	394
[226]	395	#ifdef MESOSCALE
	396	#include "meso_inc/meso_inc_caps.F"
	397	#endif
[42]	398
[226]	399	#ifndef MESOSCALE
	400	c Initialize thermospheric parameters
	401	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[42]	402
[226]	403	if (callthermos) call param_read
	404	#endif
[42]	405	c Initialize R and Cp as constant
	406
	407	if (.not.callthermos .and. .not.photochem) then
	408	do l=1,nlayermx
	409	do ig=1,ngridmx
	410	rnew(ig,l)=r
	411	cpnew(ig,l)=cpp
	412	mmean(ig,l)=mugaz
	413	enddo
	414	enddo
	415	endif
	416
	417	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
	418	write(,)"physiq: water_param Surface ice alb:",alb_surfice
	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 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	437	call zerophys(ngrid*nlayer, pdv)
	438	call zerophys(ngrid*nlayer, pdu)
	439	call zerophys(ngrid*nlayer, pdt)
	440	call zerophys(ngridnlayernq, pdq)
	441	call zerophys(ngrid, pdpsrf)
	442	call zerophys(ngrid, zflubid)
	443	call zerophys(ngrid, zdtsurf)
	444	call zerophys(ngrid*nq, dqsurf)
	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 ------------------------------------------
	550	print, 'check lat lon', lati(igout)180/pi,
	551	. long(igout)*180/pi
	552	print, 'Ls =',zls180./pi
	553	print*, 'tauref(700 Pa) =',tauref(igout)
	554	print, '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 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	598	call zerophys(ngrid*nlayer,zdtnirco2)
	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
	683	CALL zerophys(ngrid*nlayer,zdum1)
	684	CALL zerophys(ngrid*nlayer,zdum2)
	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
[42]	691	c Calling vdif (Martian version WITH CO2 condensation)
	692	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
	693	$ ptimestep,capcal,lwrite,
	694	$ pplay,pplev,zzlay,zzlev,z0,
	695	$ pu,pv,zh,pq,tsurf,emis,qsurf,
	696	$ zdum1,zdum2,zdh,pdq,zflubid,
	697	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
[256]	698	& zdqdif,zdqsdif,wmax_th,zcdv)
[42]	699
[226]	700	#ifdef MESOSCALE
	701	#include "meso_inc/meso_inc_les.F"
	702	#endif
[42]	703	DO l=1,nlayer
	704	DO ig=1,ngrid
	705	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	706	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	707	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	708
	709	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	710
	711	ENDDO
	712	ENDDO
	713
[226]	714	DO ig=1,ngrid
	715	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	716	ENDDO
[42]	717
	718	if (tracer) then
	719	DO iq=1, nq
	720	DO l=1,nlayer
	721	DO ig=1,ngrid
	722	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	723	ENDDO
	724	ENDDO
	725	ENDDO
	726	DO iq=1, nq
	727	DO ig=1,ngrid
	728	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	729	ENDDO
	730	ENDDO
	731	end if ! of if (tracer)
	732
	733	ELSE
	734	DO ig=1,ngrid
	735	zdtsurf(ig)=zdtsurf(ig)+
	736	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	737	ENDDO
[226]	738	#ifdef MESOSCALE
[42]	739	IF (flag_LES) THEN
	740	write(,) 'LES mode !'
	741	write(,) 'Please set calldifv to T in callphys.def'
	742	STOP
	743	ENDIF
[226]	744	#endif
[42]	745	ENDIF ! of IF (calldifv)
	746
[162]	747	c-----------------------------------------------------------------------
	748	c TEST. Thermals :
	749	c HIGHLY EXPERIMENTAL, BEWARE !!
	750	c -----------------------------
	751
	752	if(calltherm) then
	753
[185]	754	call calltherm_interface(firstcall,
[162]	755	$ long,lati,zzlev,zzlay,
	756	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
[185]	757	$ pplay,pplev,pphi,zpopsk,
	758	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
[226]	759	$ dtke_th,hfmax_th,wmax_th)
	760
[162]	761	DO l=1,nlayer
	762	DO ig=1,ngrid
	763	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
	764	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
	765	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
	766	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
	767	ENDDO
	768	ENDDO
	769
	770	DO ig=1,ngrid
	771	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
	772	ENDDO
	773
	774	if (tracer) then
	775	DO iq=1,nq
	776	DO l=1,nlayer
	777	DO ig=1,ngrid
	778	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
	779	ENDDO
	780	ENDDO
	781	ENDDO
	782	endif
[42]	783
[162]	784	else !of if calltherm
	785	lmax_th(:)=0
	786	end if
	787
[42]	788	c-----------------------------------------------------------------------
	789	c 5. Dry convective adjustment:
	790	c -----------------------------
	791
	792	IF(calladj) THEN
	793
	794	DO l=1,nlayer
	795	DO ig=1,ngrid
	796	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	797	ENDDO
	798	ENDDO
	799	CALL zerophys(ngrid*nlayer,zduadj)
	800	CALL zerophys(ngrid*nlayer,zdvadj)
	801	CALL zerophys(ngrid*nlayer,zdhadj)
	802
	803	CALL convadj(ngrid,nlayer,nq,ptimestep,
[162]	804	$ pplay,pplev,zpopsk,lmax_th,
[42]	805	$ pu,pv,zh,pq,
	806	$ pdu,pdv,zdh,pdq,
	807	$ zduadj,zdvadj,zdhadj,
	808	$ zdqadj)
	809
[162]	810
[42]	811	DO l=1,nlayer
	812	DO ig=1,ngrid
	813	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	814	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	815	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	816
	817	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	818	ENDDO
	819	ENDDO
	820
	821	if(tracer) then
	822	DO iq=1, nq
	823	DO l=1,nlayer
	824	DO ig=1,ngrid
	825	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	826	ENDDO
	827	ENDDO
	828	ENDDO
	829	end if
	830	ENDIF ! of IF(calladj)
	831
	832	c-----------------------------------------------------------------------
	833	c 6. Carbon dioxide condensation-sublimation:
	834	c -------------------------------------------
	835
	836	IF (callcond) THEN
[86]	837	CALL newcondens(ngrid,nlayer,nq,ptimestep,
[42]	838	$ capcal,pplay,pplev,tsurf,pt,
	839	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	840	$ co2ice,albedo,emis,
	841	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
[226]	842	$ fluxsurf_sw,zls)
[42]	843
	844	DO l=1,nlayer
	845	DO ig=1,ngrid
	846	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	847	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	848	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	849	ENDDO
	850	ENDDO
	851	DO ig=1,ngrid
	852	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	853	ENDDO
	854
	855	IF (tracer) THEN
	856	DO iq=1, nq
	857	DO l=1,nlayer
	858	DO ig=1,ngrid
	859	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	860	ENDDO
	861	ENDDO
	862	ENDDO
	863	ENDIF ! of IF (tracer)
	864
	865	ENDIF ! of IF (callcond)
	866
	867	c-----------------------------------------------------------------------
	868	c 7. Specific parameterizations for tracers
	869	c: -----------------------------------------
	870
	871	if (tracer) then
	872
	873	c 7a. Water and ice
	874	c ---------------
	875
	876	c ---------------------------------------
	877	c Water ice condensation in the atmosphere
	878	c ----------------------------------------
	879	IF (water) THEN
	880
	881	call watercloud(ngrid,nlayer,ptimestep,
	882	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
	883	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
	884	& nq,naerkind,tau,
	885	& ccn,rdust,rice,nuice)
	886	if (activice) then
	887	c Temperature variation due to latent heat release
	888	DO l=1,nlayer
	889	DO ig=1,ngrid
	890	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
	891	ENDDO
	892	ENDDO
	893	endif
	894
	895	! increment water vapour and ice atmospheric tracers tendencies
	896	IF (water) THEN
	897	DO l=1,nlayer
	898	DO ig=1,ngrid
	899	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	900	& zdqcloud(ig,l,igcm_h2o_vap)
	901	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	902	& zdqcloud(ig,l,igcm_h2o_ice)
	903	ENDDO
	904	ENDDO
	905	ENDIF ! of IF (water) THEN
	906	! Increment water ice surface tracer tendency
	907	DO ig=1,ngrid
	908	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
	909	& zdqscloud(ig,igcm_h2o_ice)
	910	ENDDO
	911
	912	END IF ! of IF (water)
	913
	914	c 7b. Chemical species
	915	c ------------------
	916
[226]	917	#ifndef MESOSCALE
	918	c --------------
	919	c photochemistry :
	920	c --------------
	921	IF (photochem .or. thermochem) then
	922	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
	923	& zzlay,zday,pq,pdq,rice,
	924	& zdqchim,zdqschim,zdqcloud,zdqscloud)
	925	!NB: Photochemistry includes condensation of H2O2
[42]	926
[226]	927	! increment values of tracers:
	928	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	929	! tracers is zero anyways
	930	DO l=1,nlayer
	931	DO ig=1,ngrid
	932	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	933	ENDDO
	934	ENDDO
	935	ENDDO ! of DO iq=1,nq
	936	! add condensation tendency for H2O2
	937	if (igcm_h2o2.ne.0) then
	938	DO l=1,nlayer
	939	DO ig=1,ngrid
	940	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	941	& +zdqcloud(ig,l,igcm_h2o2)
	942	ENDDO
	943	ENDDO
	944	endif
	945
	946	! increment surface values of tracers:
	947	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	948	! tracers is zero anyways
	949	DO ig=1,ngrid
	950	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	951	ENDDO
	952	ENDDO ! of DO iq=1,nq
	953	! add condensation tendency for H2O2
	954	if (igcm_h2o2.ne.0) then
	955	DO ig=1,ngrid
	956	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	957	& +zdqscloud(ig,igcm_h2o2)
	958	ENDDO
	959	endif
	960
	961	END IF ! of IF (photochem.or.thermochem)
	962	#endif
	963
[42]	964	c 7c. Aerosol particles
	965	c -------------------
	966
	967	c ----------
	968	c Dust devil :
	969	c ----------
	970	IF(callddevil) then
	971	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
	972	& zdqdev,zdqsdev)
	973
	974	if (dustbin.ge.1) then
	975	do iq=1,nq
	976	DO l=1,nlayer
	977	DO ig=1,ngrid
	978	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	979	ENDDO
	980	ENDDO
	981	enddo
	982	do iq=1,nq
	983	DO ig=1,ngrid
	984	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
	985	ENDDO
	986	enddo
	987	endif ! of if (dustbin.ge.1)
	988
	989	END IF ! of IF (callddevil)
	990
	991	c -------------
	992	c Sedimentation : acts also on water ice
	993	c -------------
	994	IF (sedimentation) THEN
	995	!call zerophys(ngridnlayernq, zdqsed)
	996	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	997	!call zerophys(ngrid*nq, zdqssed)
	998	zdqssed(1:ngrid,1:nq)=0
	999
	1000	call callsedim(ngrid,nlayer, ptimestep,
	1001	& pplev,zzlev, pt, rdust, rice,
	1002	& pq, pdq, zdqsed, zdqssed,nq)
	1003
	1004	DO iq=1, nq
	1005	DO l=1,nlayer
	1006	DO ig=1,ngrid
	1007	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1008	ENDDO
	1009	ENDDO
	1010	ENDDO
	1011	DO iq=1, nq
	1012	DO ig=1,ngrid
	1013	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1014	ENDDO
	1015	ENDDO
	1016	END IF ! of IF (sedimentation)
	1017
	1018	c 7d. Updates
	1019	c ---------
	1020
	1021	DO iq=1, nq
	1022	DO ig=1,ngrid
	1023
	1024	c ---------------------------------
	1025	c Updating tracer budget on surface
	1026	c ---------------------------------
	1027	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1028
	1029	ENDDO ! (ig)
	1030	ENDDO ! (iq)
	1031
	1032	endif ! of if (tracer)
	1033
[226]	1034	#ifndef MESOSCALE
	1035	c-----------------------------------------------------------------------
	1036	c 8. THERMOSPHERE CALCULATION
	1037	c-----------------------------------------------------------------------
[42]	1038
[226]	1039	if (callthermos) then
	1040	call thermosphere(pplev,pplay,dist_sol,
	1041	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	1042	& pt,pq,pu,pv,pdt,pdq,
	1043	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	1044
	1045	DO l=1,nlayer
	1046	DO ig=1,ngrid
	1047	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	1048	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	1049	& +zdteuv(ig,l)
	1050	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	1051	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	1052	DO iq=1, nq
	1053	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	1054	ENDDO
	1055	ENDDO
	1056	ENDDO
	1057
	1058	endif ! of if (callthermos)
	1059	#endif
	1060
[42]	1061	c-----------------------------------------------------------------------
	1062	c 9. Surface and sub-surface soil temperature
	1063	c-----------------------------------------------------------------------
	1064	c
	1065	c
	1066	c 9.1 Increment Surface temperature:
	1067	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1068
	1069	DO ig=1,ngrid
	1070	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1071	ENDDO
	1072
	1073	c Prescribe a cold trap at south pole (except at high obliquity !!)
	1074	c Temperature at the surface is set there to be the temperature
	1075	c corresponding to equilibrium temperature between phases of CO2
	1076
	1077	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
[226]	1078	#ifndef MESOSCALE
	1079	if (caps.and.(obliquit.lt.27.)) then
	1080	! NB: Updated surface pressure, at grid point 'ngrid', is
	1081	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
	1082	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	1083	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	1084	endif
	1085	#endif
[42]	1086	c -------------------------------------------------------------
	1087	c Change of surface albedo (set to 0.4) in case of ground frost
	1088	c everywhere except on the north permanent cap and in regions
	1089	c covered by dry ice.
	1090	c ALWAYS PLACE these lines after newcondens !!!
	1091	c -------------------------------------------------------------
	1092	do ig=1,ngrid
	1093	if ((co2ice(ig).eq.0).and.
	1094	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
	1095	albedo(ig,1) = alb_surfice
	1096	albedo(ig,2) = alb_surfice
	1097	endif
	1098	enddo ! of do ig=1,ngrid
	1099	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
	1100
	1101	c
	1102	c 9.2 Compute soil temperatures and subsurface heat flux:
	1103	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1104	IF (callsoil) THEN
	1105	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
	1106	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1107	ENDIF
	1108
	1109	c-----------------------------------------------------------------------
	1110	c 10. Write output files
	1111	c ----------------------
	1112
	1113	c -------------------------------
	1114	c Dynamical fields incrementation
	1115	c -------------------------------
	1116	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1117	! temperature, zonal and meridional wind
	1118	DO l=1,nlayer
	1119	DO ig=1,ngrid
	1120	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1121	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1122	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1123	ENDDO
	1124	ENDDO
	1125
	1126	! tracers
	1127	DO iq=1, nq
	1128	DO l=1,nlayer
	1129	DO ig=1,ngrid
	1130	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1131	ENDDO
	1132	ENDDO
	1133	ENDDO
	1134
	1135	! surface pressure
	1136	DO ig=1,ngrid
	1137	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1138	ENDDO
	1139
	1140	! pressure
	1141	DO l=1,nlayer
	1142	DO ig=1,ngrid
	1143	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1144	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1145	ENDDO
	1146	ENDDO
	1147
	1148	! Density
	1149	DO l=1,nlayer
	1150	DO ig=1,ngrid
	1151	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1152	ENDDO
	1153	ENDDO
	1154
[226]	1155	c Compute surface stress : (NB: z0 is a common in surfdat.h)
[42]	1156	c DO ig=1,ngrid
[226]	1157	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
[42]	1158	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1159	c ENDDO
	1160
	1161	c Sum of fluxes in solar spectral bands (for output only)
	1162	DO ig=1,ngrid
[226]	1163	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1164	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
[42]	1165	ENDDO
	1166	c ***** TEST ****************************************************
	1167	ztim1 = 999
	1168	DO l=1,nlayer
	1169	DO ig=1,ngrid
	1170	if (pt(ig,l).lt.ztim1) then
	1171	ztim1 = pt(ig,l)
	1172	igmin = ig
	1173	lmin = l
	1174	end if
	1175	ENDDO
	1176	ENDDO
[234]	1177	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
[42]	1178	write(,) 'PHYSIQ: stability WARNING :'
	1179	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1180	& 'ig l =', igmin, lmin
	1181	end if
	1182	c *******************************************************************
	1183
	1184	c ---------------------
	1185	c Outputs to the screen
	1186	c ---------------------
	1187
	1188	IF (lwrite) THEN
	1189	PRINT*,'Global diagnostics for the physics'
	1190	PRINT,'Variables and their increments x and dx/dt dt'
	1191	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1192	WRITE(*,'(2f10.5,2f15.5)')
	1193	s tsurf(igout),zdtsurf(igout)*ptimestep,
	1194	s pplev(igout,1),pdpsrf(igout)*ptimestep
	1195	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1196	WRITE(*,'(i4,6f10.5)') (l,
	1197	s pu(igout,l),pdu(igout,l)*ptimestep,
	1198	s pv(igout,l),pdv(igout,l)*ptimestep,
	1199	s pt(igout,l),pdt(igout,l)*ptimestep,
	1200	s l=1,nlayer)
	1201	ENDIF ! of IF (lwrite)
	1202
	1203	IF (ngrid.NE.1) THEN
	1204
[226]	1205	#ifndef MESOSCALE
[42]	1206	c -------------------------------------------------------------------
	1207	c Writing NetCDF file "RESTARTFI" at the end of the run
	1208	c -------------------------------------------------------------------
	1209	c Note: 'restartfi' is stored just before dynamics are stored
	1210	c in 'restart'. Between now and the writting of 'restart',
	1211	c there will have been the itau=itau+1 instruction and
	1212	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1213	c thus we store for time=time+dtvr
	1214
[226]	1215	IF(lastcall) THEN
	1216	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1217	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1218	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1219	. ptimestep,pday,
	1220	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
	1221	. area,albedodat,inertiedat,zmea,zstd,zsig,
	1222	. zgam,zthe)
	1223	ENDIF
	1224	#endif
[42]	1225
	1226	c -------------------------------------------------------------------
	1227	c Calculation of diagnostic variables written in both stats and
	1228	c diagfi files
	1229	c -------------------------------------------------------------------
	1230
	1231	if (tracer) then
	1232	if (water) then
	1233
	1234	call zerophys(ngrid,mtot)
	1235	call zerophys(ngrid,icetot)
	1236	call zerophys(ngrid,rave)
	1237	call zerophys(ngrid,tauTES)
	1238	do ig=1,ngrid
	1239	do l=1,nlayermx
	1240	mtot(ig) = mtot(ig) +
	1241	& zq(ig,l,igcm_h2o_vap) *
	1242	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1243	icetot(ig) = icetot(ig) +
	1244	& zq(ig,l,igcm_h2o_ice) *
	1245	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1246	rave(ig) = rave(ig) +
	1247	& zq(ig,l,igcm_h2o_ice) *
	1248	& (pplev(ig,l) - pplev(ig,l+1)) / g *
	1249	& rice(ig,l) * (1.+nuice_ref)
	1250	c Computing abs optical depth at 825 cm-1 in each
	1251	c layer to simulate NEW TES retrieval
	1252	Qabsice = min(
	1253	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1254	& )
	1255	opTES(ig,l)= 0.75 * Qabsice *
	1256	& zq(ig,l,igcm_h2o_ice) *
	1257	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1258	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1259	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1260	enddo
	1261	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
	1262	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1263	enddo
	1264
	1265	endif ! of if (water)
	1266	endif ! of if (tracer)
	1267
	1268	c -----------------------------------------------------------------
	1269	c WSTATS: Saving statistics
	1270	c -----------------------------------------------------------------
	1271	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1272	c which can later be used to make the statistic files of the run:
	1273	c "stats") only possible in 3D runs !
	1274
	1275	IF (callstats) THEN
	1276
[226]	1277	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1278	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1279	call wstats(ngrid,"co2ice","CO2 ice cover",
	1280	& "kg.m-2",2,co2ice)
	1281	call wstats(ngrid,"fluxsurf_lw",
	1282	& "Thermal IR radiative flux to surface","W.m-2",2,
	1283	& fluxsurf_lw)
	1284	call wstats(ngrid,"fluxsurf_sw",
	1285	& "Solar radiative flux to surface","W.m-2",2,
	1286	& fluxsurf_sw_tot)
	1287	call wstats(ngrid,"fluxtop_lw",
	1288	& "Thermal IR radiative flux to space","W.m-2",2,
	1289	& fluxtop_lw)
	1290	call wstats(ngrid,"fluxtop_sw",
	1291	& "Solar radiative flux to space","W.m-2",2,
	1292	& fluxtop_sw_tot)
	1293	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1294	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1295	call wstats(ngrid,"v","Meridional (North-South) wind",
	1296	& "m.s-1",3,zv)
	1297	call wstats(ngrid,"w","Vertical (down-up) wind",
	1298	& "m.s-1",3,pw)
	1299	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
	1300	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
	1301	c call wstats(ngrid,"q2",
	1302	c & "Boundary layer eddy kinetic energy",
	1303	c & "m2.s-2",3,q2)
	1304	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1305	c & emis)
	1306	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1307	c & 2,zstress)
	1308	c call wstats(ngrid,"sw_htrt","sw heat.rate",
	1309	c & "W.m-2",3,zdtsw)
	1310	c call wstats(ngrid,"lw_htrt","lw heat.rate",
	1311	c & "W.m-2",3,zdtlw)
[42]	1312
[226]	1313	if (tracer) then
	1314	if (water) then
	1315	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1316	& *mugaz/mmol(igcm_h2o_vap)
	1317	call wstats(ngrid,"vmr_h2ovapor",
	1318	& "H2O vapor volume mixing ratio","mol/mol",
	1319	& 3,vmr)
	1320	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1321	& *mugaz/mmol(igcm_h2o_ice)
	1322	call wstats(ngrid,"vmr_h2oice",
	1323	& "H2O ice volume mixing ratio","mol/mol",
	1324	& 3,vmr)
	1325	call wstats(ngrid,"h2o_ice_s",
	1326	& "surface h2o_ice","kg/m2",
	1327	& 2,qsurf(1,igcm_h2o_ice))
[42]	1328
[226]	1329	call wstats(ngrid,"mtot",
	1330	& "total mass of water vapor","kg/m2",
	1331	& 2,mtot)
	1332	call wstats(ngrid,"icetot",
	1333	& "total mass of water ice","kg/m2",
	1334	& 2,icetot)
	1335	call wstats(ngrid,"reffice",
	1336	& "Mean reff","m",
	1337	& 2,rave)
	1338	c call wstats(ngrid,"rice",
	1339	c & "Ice particle size","m",
	1340	c & 3,rice)
	1341	c If activice is true, tauTES is computed in aeropacity.F;
	1342	if (.not.activice) then
	1343	call wstats(ngrid,"tauTESap",
	1344	& "tau abs 825 cm-1","",
	1345	& 2,tauTES)
	1346	endif
	1347
	1348	endif ! of if (water)
	1349
	1350	if (thermochem.or.photochem) then
	1351	do iq=1,nq
	1352	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
	1353	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
	1354	. (noms(iq).eq."h2").or.
	1355	. (noms(iq).eq."o3")) then
	1356	do l=1,nlayer
	1357	do ig=1,ngrid
	1358	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
	1359	end do
	1360	end do
	1361	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1362	. "Volume mixing ratio","mol/mol",3,vmr)
	1363	endif
	1364	enddo
	1365	endif ! of if (thermochem.or.photochem)
	1366
	1367	endif ! of if (tracer)
	1368
	1369	IF(lastcall) THEN
	1370	write (,) "Writing stats..."
	1371	call mkstats(ierr)
	1372	ENDIF
	1373
[42]	1374	ENDIF !if callstats
	1375
	1376	c (Store EOF for Mars Climate database software)
	1377	IF (calleofdump) THEN
	1378	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1379	ENDIF
	1380
[234]	1381
[226]	1382	#ifdef MESOSCALE
[234]	1383	!!!
	1384	!!! OUTPUT FIELDS
	1385	!!!
[81]	1386	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
	1387	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
[82]	1388	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
[81]	1389	IF (igcm_h2o_ice .ne. 0) THEN
	1390	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
	1391	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	1392	. * mugaz / mmol(igcm_h2o_ice)
	1393	ENDIF
[226]	1394	c AUTOMATICALLY GENERATED FROM REGISTRY
[42]	1395	#include "fill_save.inc"
[226]	1396	#else
[42]	1397
[226]	1398	c ==========================================================
	1399	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1400	c any variable for diagnostic (output with period
	1401	c "ecritphy", set in "run.def")
	1402	c ==========================================================
	1403	c WRITEDIAGFI can ALSO be called from any other subroutines
	1404	c for any variables !!
	1405	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1406	c & emis)
	1407	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
	1408	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
	1409	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1410	& tsurf)
	1411	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1412	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1413	& co2ice)
	1414	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1415	c & "K",2,zt(1,7))
	1416	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1417	& fluxsurf_lw)
	1418	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1419	& fluxsurf_sw_tot)
	1420	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1421	& fluxtop_lw)
	1422	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1423	& fluxtop_sw_tot)
	1424	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1425	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1426	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1427	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1428	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
	1429	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
	1430	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
	1431	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
	1432	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1433	c & zstress)
	1434	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	1435	c & 'w.m-2',3,zdtsw)
	1436	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	1437	c & 'w.m-2',3,zdtlw)
	1438	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1439	c & 'tau abs 825 cm-1',
	1440	c & '',2,tauTES)
[42]	1441
	1442
[226]	1443	c ----------------------------------------------------------
	1444	c Outputs of the CO2 cycle
	1445	c ----------------------------------------------------------
[42]	1446
[226]	1447	if (tracer.and.(igcm_co2.ne.0)) then
	1448	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	1449	! & "kg/kg",2,zq(1,1,igcm_co2))
	1450	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	1451	! & "kg/kg",3,zq(1,1,igcm_co2))
	1452
	1453	! Compute co2 column
	1454	call zerophys(ngrid,co2col)
	1455	do l=1,nlayermx
	1456	do ig=1,ngrid
	1457	co2col(ig)=co2col(ig)+
	1458	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
	1459	enddo
	1460	enddo
	1461	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	1462	& co2col)
	1463	endif ! of if (tracer.and.(igcm_co2.ne.0))
[42]	1464
[226]	1465	c ----------------------------------------------------------
	1466	c Outputs of the water cycle
	1467	c ----------------------------------------------------------
	1468	if (tracer) then
	1469	if (water) then
	1470
	1471	CALL WRITEDIAGFI(ngridmx,'mtot',
	1472	& 'total mass of water vapor',
	1473	& 'kg/m2',2,mtot)
	1474	CALL WRITEDIAGFI(ngridmx,'icetot',
	1475	& 'total mass of water ice',
	1476	& 'kg/m2',2,icetot)
	1477	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1478	c & *mugaz/mmol(igcm_h2o_ice)
	1479	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
	1480	c & 'mol/mol',3,vmr)
	1481	CALL WRITEDIAGFI(ngridmx,'reffice',
	1482	& 'Mean reff',
	1483	& 'm',2,rave)
	1484	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
	1485	c & 'm',3,rice)
	1486	c If activice is true, tauTES is computed in aeropacity.F;
	1487	if (.not.activice) then
	1488	CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1489	& 'tau abs 825 cm-1',
	1490	& '',2,tauTES)
	1491	endif
	1492	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
	1493	& 'surface h2o_ice',
	1494	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
	1495	endif !(water)
	1496
	1497
	1498	if (water.and..not.photochem) then
	1499	iq=nq
	1500	c write(str2(1:2),'(i2.2)') iq
	1501	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
	1502	c & 'kg.m-2',2,zdqscloud(1,iq))
	1503	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
	1504	c & 'kg/kg',3,zdqchim(1,1,iq))
	1505	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
	1506	c & 'kg/kg',3,zdqdif(1,1,iq))
	1507	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
	1508	c & 'kg/kg',3,zdqadj(1,1,iq))
	1509	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
	1510	c & 'kg/kg',3,zdqc(1,1,iq))
	1511	endif !(water.and..not.photochem)
	1512	endif
	1513
	1514	c ----------------------------------------------------------
	1515	c Outputs of the dust cycle
	1516	c ----------------------------------------------------------
	1517
	1518	c call WRITEDIAGFI(ngridmx,'tauref',
	1519	c & 'Dust ref opt depth','NU',2,tauref)
	1520
	1521	if (tracer.and.(dustbin.ne.0)) then
	1522	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
	1523	if (doubleq) then
	1524	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
	1525	c & 'kg.m-2',2,qsurf(1,1))
	1526	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
	1527	c & 'N.m-2',2,qsurf(1,2))
	1528	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
	1529	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
	1530	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
	1531	c & 'kg.m-2.s-1',2,zdqssed(1,1))
	1532	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
	1533	& 'm',3,rdust*ref_r0)
	1534	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	1535	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
	1536	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
	1537	c & 'part/kg',3,pq(1,1,igcm_dust_number))
	1538	else
	1539	do iq=1,dustbin
	1540	write(str2(1:2),'(i2.2)') iq
	1541	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
	1542	& 'kg/kg',3,zq(1,1,iq))
	1543	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
	1544	& 'kg.m-2',2,qsurf(1,iq))
	1545	end do
	1546	endif ! (doubleq)
	1547	c if (submicron) then
	1548	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
	1549	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
	1550	c endif ! (submicron)
	1551	end if ! (tracer.and.(dustbin.ne.0))
	1552
	1553	c ----------------------------------------------------------
	1554	c Outputs of thermals
	1555	c ----------------------------------------------------------
	1556	if (calltherm) then
	1557
	1558	! call WRITEDIAGFI(ngrid,'dtke',
	1559	! & 'tendance tke thermiques','m2/s2',
	1560	! & 3,dtke_th)
	1561	! call WRITEDIAGFI(ngrid,'d_u_ajs',
	1562	! & 'tendance u thermiques','m/s',
	1563	! & 3,pdu_th*ptimestep)
	1564	! call WRITEDIAGFI(ngrid,'d_v_ajs',
	1565	! & 'tendance v thermiques','m/s',
	1566	! & 3,pdv_th*ptimestep)
	1567	! if (tracer) then
	1568	! if (nq .eq. 2) then
	1569	! call WRITEDIAGFI(ngrid,'deltaq_th',
	1570	! & 'delta q thermiques','kg/kg',
	1571	! & 3,ptimestep*pdq_th(:,:,2))
	1572	! endif
	1573	! endif
	1574
	1575	lmax_th_out(:)=real(lmax_th(:))
	1576
	1577	call WRITEDIAGFI(ngridmx,'lmax_th',
	1578	& 'hauteur du thermique','K',
	1579	& 2,lmax_th_out)
	1580	call WRITEDIAGFI(ngridmx,'hfmax_th',
	1581	& 'maximum TH heat flux','K.m/s',
	1582	& 2,hfmax_th)
	1583	call WRITEDIAGFI(ngridmx,'wmax_th',
	1584	& 'maximum TH vertical velocity','m/s',
	1585	& 2,wmax_th)
	1586
	1587	endif
	1588
	1589	c ----------------------------------------------------------
	1590	c Output in netcdf file "diagsoil.nc" for subterranean
	1591	c variables (output every "ecritphy", as for writediagfi)
	1592	c ----------------------------------------------------------
	1593
	1594	! Write soil temperature
	1595	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	1596	! & 3,tsoil)
	1597	! Write surface temperature
	1598	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	1599	! & 2,tsurf)
	1600
	1601	c ==========================================================
	1602	c END OF WRITEDIAGFI
	1603	c ==========================================================
	1604	#endif
	1605
[42]	1606	ELSE ! if(ngrid.eq.1)
	1607
	1608	print,'Ls =',zls180./pi,
	1609	& ' tauref(700 Pa) =',tauref
	1610	c ----------------------------------------------------------------------
	1611	c Output in grads file "g1d" (ONLY when using testphys1d)
	1612	c (output at every X physical timestep)
	1613	c ----------------------------------------------------------------------
	1614	c
	1615	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	1616	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	1617	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	1618
	1619	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	1620	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	1621	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	1622	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	1623
[226]	1624	! THERMALS STUFF 1D
[42]	1625
[226]	1626	if(calltherm) then
	1627
	1628	lmax_th_out(:)=real(lmax_th(:))
	1629
	1630	call WRITEDIAGFI(ngridmx,'lmax_th',
	1631	& 'hauteur du thermique','point',
	1632	& 0,lmax_th_out)
	1633	call WRITEDIAGFI(ngridmx,'hfmax_th',
	1634	& 'maximum TH heat flux','K.m/s',
	1635	& 0,hfmax_th)
	1636	call WRITEDIAGFI(ngridmx,'wmax_th',
	1637	& 'maximum TH vertical velocity','m/s',
	1638	& 0,wmax_th)
	1639
	1640
	1641	co2col(:)=0.
	1642	if (tracer) then
	1643	do l=1,nlayermx
	1644	do ig=1,ngrid
	1645	co2col(ig)=co2col(ig)+
	1646	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
	1647	enddo
	1648	enddo
	1649
	1650	end if
	1651	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
	1652	& ,'kg/m-2',0,co2col)
	1653	endif
	1654	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
	1655	& ,'m/s',1,pw)
	1656	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
	1657	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
	1658	& tsurf)
	1659
	1660	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
	1661	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
	1662	! or output in diagfi.nc (for testphys1d)
	1663	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
	1664	call WRITEDIAGFI(ngridmx,'temp','Temperature',
	1665	& 'K',1,zt)
	1666
	1667	if(tracer) then
	1668	c CALL writeg1d(ngrid,1,tau,'tau','SI')
	1669	do iq=1,nq
	1670	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	1671	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
	1672	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	1673	end do
	1674	end if
	1675
	1676	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
	1677
	1678	do l=2,nlayer-1
	1679	tmean=zt(1,l)
	1680	if(zt(1,l).ne.zt(1,l-1))
	1681	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	1682	zlocal(l)= zlocal(l-1)
	1683	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
	1684	enddo
	1685	zlocal(nlayer)= zlocal(nlayer-1)-
	1686	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
	1687	& rnew(1,nlayer)*tmean/g
	1688
[42]	1689	END IF ! if(ngrid.ne.1)
	1690
	1691	icount=icount+1
	1692	RETURN
	1693	END

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

Download in other formats:

Original Format