Context Navigation

physiq.F @ 503

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

M 498 libf/phymars/calltherm_interface.F90
^{------------------> Added computation of wstar. Wmax is no longer an output.}

M 498 libf/phymars/physiq.F
^{------------------> Minor changes.}

M 498 libf/phymars/vdif_cd.F
^{------------------> Changed computation of bulk Richardson number to allow more unstable configurations.}

M 498 libf/phymars/vdifc.F
^{------------------> Changed computation of surface conductance coefficients following a new approach. Comparisons with LES are much better.}

File size: 73.5 KB

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

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

Download in other formats:

Original Format