Context Navigation

physiq.F @ 1038

Last change on this file since 1038 was 1038, checked in by aslmd, 11 years ago

MESOSCALE
LMDZ.MARS

--> Performed the necessary modifications for dynamic tracers

to work with the mesoscale model (new physics).

--> Added precompiling flag MESOSCALE around pressure modifications

done in revision 883. This makes the mesoscale model become crazy.

--> Added an option -e in makemeso to erase a configuration and start over.

NOTE
--> not sure recent versions (rev>1000) are compliant with nesting compilation.
--> use mesoscale model + new physics with caution. still not stabilized.

File size: 92.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
[1036]	14	use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2,
	15	& igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
	16	& igcm_ccn_mass, igcm_ccn_number,
	17	& igcm_dust_mass, igcm_dust_number, igcm_h2o2,
	18	& nuice_ref, rho_ice, rho_dust, ref_r0
[544]	19
[42]	20	IMPLICIT NONE
	21	c=======================================================================
	22	c
	23	c subject:
	24	c --------
	25	c
	26	c Organisation of the physical parametrisations of the LMD
	27	c martian atmospheric general circulation model.
	28	c
	29	c The GCM can be run without or with tracer transport
	30	c depending on the value of Logical "tracer" in file "callphys.def"
	31	c Tracers may be water vapor, ice OR chemical species OR dust particles
	32	c
	33	c SEE comments in initracer.F about numbering of tracer species...
	34	c
	35	c It includes:
	36	c
	37	c 1. Initialization:
	38	c 1.1 First call initializations
	39	c 1.2 Initialization for every call to physiq
	40	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	41	c 2. Compute radiative transfer tendencies
	42	c (longwave and shortwave) for CO2 and aerosols.
	43	c 3. Gravity wave and subgrid scale topography drag :
	44	c 4. Vertical diffusion (turbulent mixing):
	45	c 5. Convective adjustment
	46	c 6. Condensation and sublimation of carbon dioxide.
	47	c 7. TRACERS :
	48	c 7a. water and water ice
	49	c 7b. call for photochemistry when tracers are chemical species
	50	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
	51	c 7d. updates (CO2 pressure variations, surface budget)
	52	c 8. Contribution to tendencies due to thermosphere
	53	c 9. Surface and sub-surface temperature calculations
	54	c 10. Write outputs :
	55	c - "startfi", "histfi" (if it's time)
	56	c - Saving statistics (if "callstats = .true.")
	57	c - Dumping eof (if "calleofdump = .true.")
	58	c - Output any needed variables in "diagfi"
	59	c 11. Diagnostic: mass conservation of tracers
	60	c
	61	c author:
	62	c -------
	63	c Frederic Hourdin 15/10/93
	64	c Francois Forget 1994
	65	c Christophe Hourdin 02/1997
	66	c Subroutine completly rewritten by F.Forget (01/2000)
	67	c Introduction of the photochemical module: S. Lebonnois (11/2002)
	68	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
	69	c Water ice clouds: Franck Montmessin (update 06/2003)
	70	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	71	c Nb: See callradite.F for more information.
[234]	72	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
[414]	73	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
[42]	74	c
	75	c arguments:
	76	c ----------
	77	c
	78	c input:
	79	c ------
	80	c ecri period (in dynamical timestep) to write output
	81	c ngrid Size of the horizontal grid.
	82	c All internal loops are performed on that grid.
	83	c nlayer Number of vertical layers.
	84	c nq Number of advected fields
	85	c firstcall True at the first call
	86	c lastcall True at the last call
	87	c pday Number of days counted from the North. Spring
	88	c equinoxe.
	89	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	90	c ptimestep timestep (s)
	91	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	92	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	93	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	94	c pu(ngrid,nlayer) u component of the wind (ms-1)
	95	c pv(ngrid,nlayer) v component of the wind (ms-1)
	96	c pt(ngrid,nlayer) Temperature (K)
	97	c pq(ngrid,nlayer,nq) Advected fields
[330]	98	c pudyn(ngrid,nlayer) \|
	99	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
	100	c ptdyn(ngrid,nlayer) \| corresponding variables
	101	c pqdyn(ngrid,nlayer,nq) \|
[42]	102	c pw(ngrid,?) vertical velocity
	103	c
	104	c output:
	105	c -------
	106	c
[330]	107	c pdu(ngrid,nlayermx) \|
	108	c pdv(ngrid,nlayermx) \| Temporal derivative of the corresponding
	109	c pdt(ngrid,nlayermx) \| variables due to physical processes.
[1036]	110	c pdq(ngrid,nlayermx,nq) \|
[330]	111	c pdpsrf(ngrid) \|
[42]	112	c tracerdyn call tracer in dynamical part of GCM ?
	113
	114	c
	115	c=======================================================================
	116	c
	117	c 0. Declarations :
	118	c ------------------
	119
	120	#include "dimensions.h"
	121	#include "dimphys.h"
	122	#include "comgeomfi.h"
	123	#include "surfdat.h"
[226]	124	#include "comsoil.h"
[42]	125	#include "comdiurn.h"
	126	#include "callkeys.h"
	127	#include "comcstfi.h"
	128	#include "planete.h"
	129	#include "comsaison.h"
	130	#include "control.h"
	131	#include "dimradmars.h"
	132	#include "comg1d.h"
[1036]	133	!#include "tracer.h"
[42]	134	#include "nlteparams.h"
[883]	135	#include "comvert.h"
[42]	136
	137	#include "chimiedata.h"
	138	#include "param.h"
[635]	139	#include "param_v4.h"
[42]	140	#include "conc.h"
	141
	142	#include "netcdf.inc"
	143
[234]	144	#include "slope.h"
	145
[86]	146	#ifdef MESOSCALE
[42]	147	#include "wrf_output_2d.h"
	148	#include "wrf_output_3d.h"
[1036]	149	!#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
[226]	150	#include "meso_inc/meso_inc_var.F"
[86]	151	#endif
[42]	152
	153	c Arguments :
	154	c -----------
	155
	156	c inputs:
	157	c -------
	158	INTEGER ngrid,nlayer,nq
	159	REAL ptimestep
	160	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	161	REAL pphi(ngridmx,nlayer)
	162	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	163	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	164	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
	165	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	166	LOGICAL firstcall,lastcall
[226]	167
[42]	168	REAL pday
[999]	169	REAL ptime
[42]	170	logical tracerdyn
	171
	172	c outputs:
	173	c --------
	174	c physical tendencies
	175	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	176	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	177	REAL pdpsrf(ngridmx) ! surface pressure tendency
	178
	179
	180	c Local saved variables:
	181	c ----------------------
	182	c aerosol (dust or ice) extinction optical depth at reference wavelength
	183	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
	184	c aerosol optical properties :
	185	REAL aerosol(ngridmx,nlayermx,naerkind)
	186
	187	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	188	INTEGER icount ! counter of calls to physiq during the run.
	189	REAL tsurf(ngridmx) ! Surface temperature (K)
	190	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	191	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
	192	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
	193	REAL emis(ngridmx) ! Thermal IR surface emissivity
	194	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	195	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	196	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	197	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	198	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
[1036]	199	REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
[283]	200	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	201
[42]	202	c Variables used by the water ice microphysical scheme:
	203	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
	204	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
	205	! of the size distribution
[358]	206	real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
	207	real rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
[576]	208	REAL surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3, if photochemistry)
	209	REAL surfice(ngridmx,nlayermx) ! ice surface area (m2/m3, if photochemistry)
[833]	210	REAL inertiesoil(ngridmx,nsoilmx)! Time varying subsurface
	211	! thermal inertia (J.s-1/2.m-2.K-1)
	212	! (used only when tifeedback=.true.)
[42]	213
[234]	214	c Variables used by the slope model
	215	REAL sl_ls, sl_lct, sl_lat
	216	REAL sl_tau, sl_alb, sl_the, sl_psi
	217	REAL sl_fl0, sl_flu
	218	REAL sl_ra, sl_di0
	219	REAL sky
	220
[999]	221	SAVE day_ini, icount, time_phys
[42]	222	SAVE aerosol, tsurf,tsoil
	223	SAVE co2ice,albedo,emis, q2
[1036]	224	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky
[42]	225
	226	REAL stephan
	227	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	228	SAVE stephan
	229
	230	c Local variables :
	231	c -----------------
	232
	233	REAL CBRT
	234	EXTERNAL CBRT
	235
[1036]	236	! CHARACTER*80 fichier
	237	INTEGER l,ig,ierr,igout,iq,tapphys
[42]	238
	239	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
	240	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
	241	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
	242	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
[627]	243	REAL tauref(ngridmx) ! Reference column optical depth at odpref
	244	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
[42]	245	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	246	REAL zls ! solar longitude (rad)
	247	REAL zday ! date (time since Ls=0, in martian days)
	248	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	249	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
[1036]	250	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
[42]	251
	252	c Tendancies due to various processes:
[1036]	253	REAL dqsurf(ngridmx,nq)
[42]	254	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	255	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
[1036]	256	! REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
	257	! REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
[42]	258	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
	259	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
	260	REAL zdtsurf(ngridmx) ! (K/s)
	261	REAL zdtcloud(ngridmx,nlayermx)
	262	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	263	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	264	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	265	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	266	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	267	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	268	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	269	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	270
[1036]	271	REAL zdqdif(ngridmx,nlayermx,nq), zdqsdif(ngridmx,nq)
	272	REAL zdqsed(ngridmx,nlayermx,nq), zdqssed(ngridmx,nq)
	273	REAL zdqdev(ngridmx,nlayermx,nq), zdqsdev(ngridmx,nq)
	274	REAL zdqadj(ngridmx,nlayermx,nq)
	275	REAL zdqc(ngridmx,nlayermx,nq)
	276	REAL zdqcloud(ngridmx,nlayermx,nq)
	277	REAL zdqscloud(ngridmx,nq)
	278	REAL zdqchim(ngridmx,nlayermx,nq)
	279	REAL zdqschim(ngridmx,nq)
[42]	280
	281	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	282	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	283	REAL zdumolvis(ngridmx,nlayermx)
	284	REAL zdvmolvis(ngridmx,nlayermx)
[1036]	285	real zdqmoldiff(ngridmx,nlayermx,nq)
[42]	286
	287	c Local variable for local intermediate calcul:
	288	REAL zflubid(ngridmx)
	289	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	290	REAL zdum1(ngridmx,nlayermx)
	291	REAL zdum2(ngridmx,nlayermx)
	292	REAL ztim1,ztim2,ztim3, z1,z2
	293	REAL ztime_fin
	294	REAL zdh(ngridmx,nlayermx)
	295	INTEGER length
	296	PARAMETER (length=100)
	297
	298	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	299	c -----------------------------------------------------------------------------
	300	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	301	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
[1036]	302	REAL zq(ngridmx,nlayermx,nq)
[42]	303	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	304	character*2 str2
[1036]	305	! character*5 str5
[42]	306	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
[358]	307	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
	308	SAVE tauscaling ! in case iradia NE 1
[42]	309	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
	310	integer igmin, lmin
	311	logical tdiag
	312
	313	real co2col(ngridmx) ! CO2 column
[883]	314	! pplev and pplay are dynamical inputs and must not be modified in the physics.
	315	! instead, use zplay and zplev :
	316	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
[1036]	317	! REAL zstress(ngrid),cd
	318	real tmean, zlocal(nlayermx)
[42]	319	real rho(ngridmx,nlayermx) ! density
	320	real vmr(ngridmx,nlayermx) ! volume mixing ratio
[705]	321	real rhopart(ngridmx,nlayermx) ! number density of a given species
[1036]	322	real colden(ngridmx,nq) ! vertical column of tracers
[42]	323	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
	324	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
[669]	325	REAL Nccntot(ngridmx) ! Total number of ccn (nbr/m2)
	326	REAL Mccntot(ngridmx) ! Total mass of ccn (kg/m2)
[42]	327	REAL rave(ngridmx) ! Mean water ice effective radius (m)
	328	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
	329	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
	330	REAL Qabsice ! Water ice absorption coefficient
[520]	331	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
	332	! reference wavelength using
	333	! Qabs instead of Qext
	334	! (direct comparison with TES)
[719]	335
	336	REAL dqdustsurf(ngridmx) ! surface q dust flux (kg/m2/s)
	337	REAL dndustsurf(ngridmx) ! surface n dust flux (number/m2/s)
[756]	338	REAL ndust(ngridmx,nlayermx) ! true n dust (kg/kg)
	339	REAL qdust(ngridmx,nlayermx) ! true q dust (kg/kg)
	340	REAL nccn(ngridmx,nlayermx) ! true n ccn (kg/kg)
	341	REAL qccn(ngridmx,nlayermx) ! true q ccn (kg/kg)
[42]	342
[411]	343	c Test 1d/3d scavenging
[520]	344	real h2otot(ngridmx)
[719]	345	REAL satu(ngridmx,nlayermx) ! satu ratio for output
	346	REAL zqsat(ngridmx,nlayermx) ! saturation
[42]	347
	348	REAL time_phys
	349
[414]	350	! Added for new NLTE scheme
	351
	352	real co2vmr_gcm(ngridmx,nlayermx)
	353	real n2vmr_gcm(ngridmx,nlayermx)
	354	real ovmr_gcm(ngridmx,nlayermx)
	355	real covmr_gcm(ngridmx,nlayermx)
	356
	357
[267]	358	c Variables for PBL
[544]	359	REAL zz1(ngridmx)
[185]	360	REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
[528]	361	REAL, SAVE :: wstar(ngridmx)
	362	REAL, SAVE :: hfmax_th(ngridmx)
[185]	363	REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
[1036]	364	REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nq)
[566]	365	INTEGER lmax_th(ngridmx),dimout,n_out,n
	366	CHARACTER(50) zstring
[185]	367	REAL dtke_th(ngridmx,nlayermx+1)
[268]	368	REAL zcdv(ngridmx), zcdh(ngridmx)
[636]	369	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
[566]	370	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
[1038]	371	REAL u_out1(ngridmx)
[1036]	372	REAL T_out1(ngridmx)
[566]	373	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
[267]	374	REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp
[657]	375	REAL L_mo(ngridmx),vhf(ngridmx),vvv(ngridmx)
[1036]	376	! REAL zu2(ngridmx)
	377	REAL sensibFlux(ngridmx)
[790]	378
[42]	379	c=======================================================================
	380
	381	c 1. Initialisation:
	382	c -----------------
	383
	384	c 1.1 Initialisation only at first call
	385	c ---------------------------------------
	386	IF (firstcall) THEN
	387
[1036]	388	! allocate local (saved) arrays:
	389	allocate(qsurf(ngrid,nq))
	390
[42]	391	c variables set to 0
	392	c ~~~~~~~~~~~~~~~~~~
[286]	393	aerosol(:,:,:)=0
	394	dtrad(:,:)=0
[674]	395
	396	#ifndef MESOSCALE
[286]	397	fluxrad(:)=0
[528]	398	wstar(:)=0.
[674]	399	#else
	400	#include "meso_inc/meso_inc_ini_restart.F"
	401	#endif
[268]	402
[42]	403	c read startfi
	404	c ~~~~~~~~~~~~
[226]	405	#ifndef MESOSCALE
	406	! Read netcdf initial physical parameters.
	407	CALL phyetat0 ("startfi.nc",0,0,
	408	& nsoilmx,nq,
	409	& day_ini,time_phys,
	410	& tsurf,tsoil,emis,q2,qsurf,co2ice)
	411	#else
	412	#include "meso_inc/meso_inc_ini.F"
	413	#endif
[185]	414
[42]	415	if (pday.ne.day_ini) then
	416	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	417	& "physics and dynamics"
	418	write(,) "dynamics day: ",pday
	419	write(,) "physics day: ",day_ini
	420	stop
	421	endif
	422
	423	write (,) 'In physiq day_ini =', day_ini
	424
[286]	425	c initialize tracers
	426	c ~~~~~~~~~~~~~~~~~~
	427	tracerdyn=tracer
	428	IF (tracer) THEN
[1036]	429	CALL initracer(ngrid,nq,qsurf,co2ice)
[286]	430	ENDIF ! end tracer
	431
[42]	432	c Initialize albedo and orbital calculation
	433	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	434	CALL surfini(ngrid,co2ice,qsurf,albedo)
	435	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	436
	437	c initialize soil
	438	c ~~~~~~~~~~~~~~~
	439	IF (callsoil) THEN
[833]	440	c Thermal inertia feedback:
	441	IF (tifeedback) THEN
	442	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
	443	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
	444	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	445	ELSE
	446	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	447	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	448	ENDIF ! of IF (tifeedback)
[42]	449	ELSE
	450	PRINT*,
	451	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	452	DO ig=1,ngrid
	453	capcal(ig)=1.e5
	454	fluxgrd(ig)=0.
	455	ENDDO
	456	ENDIF
	457	icount=1
	458
[226]	459	#ifndef MESOSCALE
	460	c Initialize thermospheric parameters
	461	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[42]	462
[705]	463	if (callthermos) then
	464	if(solvarmod.eq.0) call param_read
	465	if(solvarmod.eq.1) call param_read_e107
	466	endif
[226]	467	#endif
[42]	468	c Initialize R and Cp as constant
	469
	470	if (.not.callthermos .and. .not.photochem) then
	471	do l=1,nlayermx
	472	do ig=1,ngridmx
	473	rnew(ig,l)=r
	474	cpnew(ig,l)=cpp
	475	mmean(ig,l)=mugaz
	476	enddo
	477	enddo
	478	endif
	479
[757]	480	if(callnlte.and.nltemodel.eq.2) call nlte_setup
[414]	481	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
[1013]	482	if(thermochem) call chemthermos_readini
[414]	483
[42]	484	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
[283]	485	write(,)"physiq: water_param Surface water ice albedo:",
	486	. albedo_h2o_ice
[42]	487	ENDIF
[900]	488
	489	#ifndef MESOSCALE
	490	if (callslope) call getslopes(phisfi)
[999]	491
	492	call physdem0("restartfi.nc",long,lati,nsoilmx,nq,
	493	. ptimestep,pday,time_phys,area,
	494	. albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
	495
[900]	496	#endif
[42]	497
	498	ENDIF ! (end of "if firstcall")
	499
[414]	500
[42]	501	c ---------------------------------------------------
	502	c 1.2 Initializations done at every physical timestep:
	503	c ---------------------------------------------------
	504	c
	505	IF (ngrid.NE.ngridmx) THEN
	506	PRINT*,'STOP in PHYSIQ'
	507	PRINT*,'Probleme de dimensions :'
	508	PRINT*,'ngrid = ',ngrid
	509	PRINT*,'ngridmx = ',ngridmx
	510	STOP
	511	ENDIF
	512
	513	c Initialize various variables
	514	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[286]	515	pdv(:,:)=0
	516	pdu(:,:)=0
	517	pdt(:,:)=0
	518	pdq(:,:,:)=0
	519	pdpsrf(:)=0
	520	zflubid(:)=0
	521	zdtsurf(:)=0
	522	dqsurf(:,:)=0
[42]	523	igout=ngrid/2+1
	524
	525
	526	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	527
	528	c Compute Solar Longitude (Ls) :
	529	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	530	if (season) then
	531	call solarlong(zday,zls)
	532	else
	533	call solarlong(float(day_ini),zls)
	534	end if
	535
[883]	536	c Initialize pressure levels
	537	c ~~~~~~~~~~~~~~~~~~
	538	zplev(:,:) = pplev(:,:)
	539	zplay(:,:) = pplay(:,:)
	540	ps(:) = pplev(:,1)
	541
	542
[42]	543	c Compute geopotential at interlayers
	544	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	545	c ponderation des altitudes au niveau des couches en dp/p
	546
	547	DO l=1,nlayer
	548	DO ig=1,ngrid
	549	zzlay(ig,l)=pphi(ig,l)/g
	550	ENDDO
	551	ENDDO
	552	DO ig=1,ngrid
	553	zzlev(ig,1)=0.
	554	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	555	ENDDO
	556	DO l=2,nlayer
	557	DO ig=1,ngrid
[883]	558	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
	559	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
[42]	560	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	561	ENDDO
	562	ENDDO
	563
	564
	565	! Potential temperature calculation not the same in physiq and dynamic
	566
	567	c Compute potential temperature
	568	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	569	DO l=1,nlayer
	570	DO ig=1,ngrid
[883]	571	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
[42]	572	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	573	ENDDO
	574	ENDDO
	575
[226]	576	#ifndef MESOSCALE
	577	c-----------------------------------------------------------------------
	578	c 1.2.5 Compute mean mass, cp, and R
	579	c --------------------------------
[42]	580
[226]	581	if(photochem.or.callthermos) then
[1036]	582	call concentrations(nq,zplay,pt,pdt,pq,pdq,ptimestep)
[226]	583	endif
	584	#endif
[42]	585	c-----------------------------------------------------------------------
	586	c 2. Compute radiative tendencies :
	587	c------------------------------------
	588
	589
	590	IF (callrad) THEN
	591	IF( MOD(icount-1,iradia).EQ.0) THEN
	592
	593	c Local Solar zenith angle
	594	c ~~~~~~~~~~~~~~~~~~~~~~~~
	595	CALL orbite(zls,dist_sol,declin)
	596
	597	IF(diurnal) THEN
	598	ztim1=SIN(declin)
	599	ztim2=COS(declin)COS(2.pi*(zday-.5))
	600	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	601
	602	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	603	s ztim1,ztim2,ztim3, mu0,fract)
	604
	605	ELSE
	606	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
	607	ENDIF
	608
	609	c NLTE cooling from CO2 emission
	610	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[414]	611	IF(callnlte) then
	612	if(nltemodel.eq.0.or.nltemodel.eq.1) then
[883]	613	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
[414]	614	else if(nltemodel.eq.2) then
[705]	615	co2vmr_gcm(1:ngrid,1:nlayer)=
	616	& pq(1:ngrid,1:nlayer,igcm_co2)*
	617	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
	618	n2vmr_gcm(1:ngrid,1:nlayer)=
	619	& pq(1:ngrid,1:nlayer,igcm_n2)*
	620	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
	621	covmr_gcm(1:ngrid,1:nlayer)=
	622	& pq(1:ngrid,1:nlayer,igcm_co)*
	623	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
	624	ovmr_gcm(1:ngrid,1:nlayer)=
	625	& pq(1:ngrid,1:nlayer,igcm_o)*
	626	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
[414]	627
[883]	628	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
[414]	629	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
	630	$ ovmr_gcm, zdtnlte )
[42]	631
[705]	632	zdtnlte(1:ngrid,1:nlayer)=
	633	& zdtnlte(1:ngrid,1:nlayer)/86400.
[414]	634	endif
[528]	635	else
	636	zdtnlte(:,:)=0.
[414]	637	endif
[42]	638
	639	c Find number of layers for LTE radiation calculations
	640	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
[883]	641	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
[42]	642
	643	c Note: Dustopacity.F has been transferred to callradite.F
	644
	645	c Call main radiative transfer scheme
	646	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	647	c Transfer through CO2 (except NIR CO2 absorption)
	648	c and aerosols (dust and water ice)
	649
	650	c Radiative transfer
	651	c ------------------
[833]	652
[42]	653	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
[883]	654	$ emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
[42]	655	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
[520]	656	$ tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
	657	$ nuice,co2ice)
[42]	658
[234]	659	c Outputs for basic check (middle of domain)
	660	c ------------------------------------------
[627]	661	write(*,'("Ls =",f11.6," check lat =",f10.6,
	662	& " lon =",f11.6)')
	663	& zls180./pi,lati(igout)180/pi,long(igout)*180/pi
	664	write(*,'(" tauref(",f4.0," Pa) =",f9.6,
	665	& " tau(",f4.0," Pa) =",f9.6)')
	666	& odpref,tauref(igout),
[883]	667	& odpref,tau(igout,1)*odpref/zplev(igout,1)
[234]	668	c ---------------------------------------------------------
	669	c Call slope parameterization for direct and scattered flux
	670	c ---------------------------------------------------------
	671	IF(callslope) THEN
	672	print *, 'Slope scheme is on and computing...'
	673	DO ig=1,ngrid
	674	sl_the = theta_sl(ig)
	675	IF (sl_the .ne. 0.) THEN
	676	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
	677	DO l=1,2
	678	sl_lct = ptime24. + 180.long(ig)/pi/15.
	679	sl_ra = pi*(1.0-sl_lct/12.)
	680	sl_lat = 180.*lati(ig)/pi
[577]	681	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
[234]	682	sl_alb = albedo(ig,l)
	683	sl_psi = psi_sl(ig)
	684	sl_fl0 = fluxsurf_sw(ig,l)
	685	sl_di0 = 0.
	686	if (mu0(ig) .gt. 0.) then
	687	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
	688	sl_di0 = sl_di0*1370./dist_sol/dist_sol
	689	sl_di0 = sl_di0/ztim1
	690	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
	691	endif
	692	! you never know (roundup concern...)
	693	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
	694	!!!!!!!!!!!!!!!!!!!!!!!!!!
	695	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
	696	& sl_tau, sl_alb, sl_the, sl_psi,
	697	& sl_di0, sl_fl0, sl_flu )
	698	!!!!!!!!!!!!!!!!!!!!!!!!!!
	699	fluxsurf_sw(ig,l) = sl_flu
	700	ENDDO
	701	!!! compute correction on IR flux as well
	702	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
	703	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
	704	ENDIF
	705	ENDDO
	706	ENDIF
	707
[42]	708	c CO2 near infrared absorption
	709	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[286]	710	zdtnirco2(:,:)=0
[42]	711	if (callnirco2) then
[883]	712	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
[42]	713	. mu0,fract,declin, zdtnirco2)
	714	endif
	715
	716	c Radiative flux from the sky absorbed by the surface (W.m-2)
	717	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	718	DO ig=1,ngrid
	719	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	720	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
	721	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
	722	ENDDO
	723
	724
	725	c Net atmospheric radiative heating rate (K.s-1)
	726	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	727	IF(callnlte) THEN
[883]	728	CALL blendrad(ngrid, nlayer, zplay,
[42]	729	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
	730	ELSE
	731	DO l=1,nlayer
	732	DO ig=1,ngrid
	733	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	734	& +zdtnirco2(ig,l)
	735	ENDDO
	736	ENDDO
	737	ENDIF
	738
	739	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	740
	741	c Transformation of the radiative tendencies:
	742	c -------------------------------------------
	743
	744	c Net radiative surface flux (W.m-2)
	745	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	746	c
	747	DO ig=1,ngrid
	748	zplanck(ig)=tsurf(ig)*tsurf(ig)
	749	zplanck(ig)=emis(ig)*
	750	$ stephanzplanck(ig)zplanck(ig)
	751	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
[234]	752	IF(callslope) THEN
	753	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
	754	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
	755	ENDIF
[42]	756	ENDDO
	757
	758	DO l=1,nlayer
	759	DO ig=1,ngrid
	760	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	761	ENDDO
	762	ENDDO
	763
	764	ENDIF ! of IF (callrad)
	765
[226]	766	c-----------------------------------------------------------------------
	767	c 3. Gravity wave and subgrid scale topography drag :
	768	c -------------------------------------------------
[42]	769
[226]	770
	771	IF(calllott)THEN
	772
	773	CALL calldrag_noro(ngrid,nlayer,ptimestep,
[883]	774	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
[226]	775
	776	DO l=1,nlayer
	777	DO ig=1,ngrid
	778	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
	779	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
	780	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
	781	ENDDO
	782	ENDDO
	783	ENDIF
[234]	784
[42]	785	c-----------------------------------------------------------------------
	786	c 4. Vertical diffusion (turbulent mixing):
	787	c -----------------------------------------
[226]	788
[42]	789	IF (calldifv) THEN
	790
	791	DO ig=1,ngrid
	792	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	793	ENDDO
	794
[286]	795	zdum1(:,:)=0
	796	zdum2(:,:)=0
[42]	797	do l=1,nlayer
	798	do ig=1,ngrid
	799	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	800	enddo
	801	enddo
[226]	802
[288]	803
	804	#ifdef MESOSCALE
[289]	805	IF (.not.flag_LES) THEN
[288]	806	#endif
	807	c ----------------------
[284]	808	c Treatment of a special case : using new surface layer (Richardson based)
	809	c without using the thermals in gcm and mesoscale can yield problems in
	810	c weakly unstable situations when winds are near to 0. For those cases, we add
	811	c a unit subgrid gustiness. Remember that thermals should be used we using the
	812	c Richardson based surface layer model.
[288]	813	IF ( .not.calltherm .and. callrichsl ) THEN
[284]	814	DO ig=1, ngridmx
	815	IF (zh(ig,1) .lt. tsurf(ig)) THEN
[528]	816	wstar(ig)=1.
	817	hfmax_th(ig)=0.2
	818	ELSE
	819	wstar(ig)=0.
	820	hfmax_th(ig)=0.
	821	ENDIF
[284]	822	ENDDO
	823	ENDIF
[288]	824	c ----------------------
	825	#ifdef MESOSCALE
[284]	826	ENDIF
	827	#endif
	828
[544]	829	IF (tke_heat_flux .ne. 0.) THEN
	830	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
[883]	831	& (-alog(zplay(:,1)/zplev(:,1)))
[544]	832	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
	833	ENDIF
[288]	834
[42]	835	c Calling vdif (Martian version WITH CO2 condensation)
	836	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
	837	$ ptimestep,capcal,lwrite,
[883]	838	$ zplay,zplev,zzlay,zzlev,z0,
[42]	839	$ pu,pv,zh,pq,tsurf,emis,qsurf,
	840	$ zdum1,zdum2,zdh,pdq,zflubid,
	841	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
[660]	842	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux
[529]	843	#ifdef MESOSCALE
	844	& ,flag_LES
	845	#endif
	846	& )
[42]	847
[529]	848
[226]	849	#ifdef MESOSCALE
	850	#include "meso_inc/meso_inc_les.F"
	851	#endif
[42]	852	DO l=1,nlayer
	853	DO ig=1,ngrid
	854	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	855	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	856	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	857
	858	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	859
	860	ENDDO
	861	ENDDO
	862
[226]	863	DO ig=1,ngrid
	864	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	865	ENDDO
[42]	866
[790]	867	if (tracer) then
	868	#ifdef MESOSCALE
	869	IF (.not.flag_LES) THEN
	870	#endif
[42]	871	DO iq=1, nq
	872	DO l=1,nlayer
	873	DO ig=1,ngrid
	874	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	875	ENDDO
	876	ENDDO
	877	ENDDO
	878	DO iq=1, nq
	879	DO ig=1,ngrid
	880	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	881	ENDDO
	882	ENDDO
[790]	883	#ifdef MESOSCALE
	884	ENDIF
	885	#endif
[42]	886	end if ! of if (tracer)
	887
	888	ELSE
	889	DO ig=1,ngrid
	890	zdtsurf(ig)=zdtsurf(ig)+
	891	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	892	ENDDO
[226]	893	#ifdef MESOSCALE
[42]	894	IF (flag_LES) THEN
	895	write(,) 'LES mode !'
	896	write(,) 'Please set calldifv to T in callphys.def'
	897	STOP
	898	ENDIF
[226]	899	#endif
[42]	900	ENDIF ! of IF (calldifv)
	901
[162]	902	c-----------------------------------------------------------------------
[566]	903	c 5. Thermals :
[162]	904	c -----------------------------
[566]	905
[162]	906	if(calltherm) then
	907
[1032]	908	call calltherm_interface(ngrid,nlayer,nq,
	909	$ tracer,igcm_co2,
[652]	910	$ zzlev,zzlay,
[162]	911	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
[883]	912	$ zplay,zplev,pphi,zpopsk,
[185]	913	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
[660]	914	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
[633]	915
[162]	916	DO l=1,nlayer
	917	DO ig=1,ngrid
	918	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
	919	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
	920	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
	921	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
	922	ENDDO
	923	ENDDO
	924
	925	DO ig=1,ngrid
	926	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
	927	ENDDO
	928
	929	if (tracer) then
	930	DO iq=1,nq
	931	DO l=1,nlayer
	932	DO ig=1,ngrid
	933	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
	934	ENDDO
	935	ENDDO
	936	ENDDO
	937	endif
[42]	938
[277]	939	lmax_th_out(:)=real(lmax_th(:))
	940
[162]	941	else !of if calltherm
	942	lmax_th(:)=0
[528]	943	wstar(:)=0.
	944	hfmax_th(:)=0.
[277]	945	lmax_th_out(:)=0.
[162]	946	end if
	947
[42]	948	c-----------------------------------------------------------------------
	949	c 5. Dry convective adjustment:
	950	c -----------------------------
	951
	952	IF(calladj) THEN
	953
	954	DO l=1,nlayer
	955	DO ig=1,ngrid
	956	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	957	ENDDO
	958	ENDDO
[286]	959	zduadj(:,:)=0
	960	zdvadj(:,:)=0
	961	zdhadj(:,:)=0
[42]	962
	963	CALL convadj(ngrid,nlayer,nq,ptimestep,
[883]	964	$ zplay,zplev,zpopsk,lmax_th,
[42]	965	$ pu,pv,zh,pq,
	966	$ pdu,pdv,zdh,pdq,
	967	$ zduadj,zdvadj,zdhadj,
	968	$ zdqadj)
	969
[162]	970
[42]	971	DO l=1,nlayer
	972	DO ig=1,ngrid
	973	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	974	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	975	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	976
	977	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	978	ENDDO
	979	ENDDO
	980
	981	if(tracer) then
	982	DO iq=1, nq
	983	DO l=1,nlayer
	984	DO ig=1,ngrid
	985	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	986	ENDDO
	987	ENDDO
	988	ENDDO
	989	end if
	990	ENDIF ! of IF(calladj)
	991
	992	c-----------------------------------------------------------------------
	993	c 6. Carbon dioxide condensation-sublimation:
	994	c -------------------------------------------
	995
[528]	996	IF (tituscap) THEN
	997	!!! get the actual co2 seasonal cap from Titus observations
	998	CALL geticecover( ngrid, 180.*zls/pi,
[333]	999	. 180.long/pi, 180.lati/pi, co2ice )
[528]	1000	co2ice = co2ice * 10000.
	1001	ENDIF
[883]	1002
	1003
	1004	pdpsrf(:) = 0
[332]	1005
[42]	1006	IF (callcond) THEN
[86]	1007	CALL newcondens(ngrid,nlayer,nq,ptimestep,
[883]	1008	$ capcal,zplay,zplev,tsurf,pt,
[42]	1009	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	1010	$ co2ice,albedo,emis,
	1011	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
[226]	1012	$ fluxsurf_sw,zls)
[42]	1013
	1014	DO l=1,nlayer
	1015	DO ig=1,ngrid
	1016	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	1017	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	1018	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	1019	ENDDO
	1020	ENDDO
	1021	DO ig=1,ngrid
	1022	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	1023	ENDDO
	1024
	1025	IF (tracer) THEN
	1026	DO iq=1, nq
	1027	DO l=1,nlayer
	1028	DO ig=1,ngrid
	1029	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	1030	ENDDO
	1031	ENDDO
	1032	ENDDO
	1033	ENDIF ! of IF (tracer)
	1034
[1038]	1035	#ifndef MESOSCALE
[883]	1036	! update surface pressure
	1037	DO ig=1,ngrid
	1038	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
	1039	ENDDO
	1040
	1041	! update pressure levels
	1042	DO l=1,nlayer
	1043	DO ig=1,ngrid
	1044	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
	1045	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
	1046	ENDDO
	1047	ENDDO
[885]	1048	zplev(:,nlayer+1) = 0.
[883]	1049
	1050	! update layers altitude
	1051	DO l=2,nlayer
	1052	DO ig=1,ngrid
	1053	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
	1054	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
	1055	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	1056	ENDDO
	1057	ENDDO
[1038]	1058	#endif
[883]	1059
[42]	1060	ENDIF ! of IF (callcond)
[883]	1061
[42]	1062
[883]	1063
[42]	1064	c-----------------------------------------------------------------------
	1065	c 7. Specific parameterizations for tracers
	1066	c: -----------------------------------------
	1067
	1068	if (tracer) then
	1069
	1070	c 7a. Water and ice
	1071	c ---------------
	1072
	1073	c ---------------------------------------
	1074	c Water ice condensation in the atmosphere
	1075	c ----------------------------------------
	1076	IF (water) THEN
	1077
	1078	call watercloud(ngrid,nlayer,ptimestep,
[883]	1079	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
[626]	1080	& pq,pdq,zdqcloud,zdtcloud,
[358]	1081	& nq,tau,tauscaling,rdust,rice,nuice,
	1082	& rsedcloud,rhocloud)
[833]	1083
[633]	1084	c Temperature variation due to latent heat release
[42]	1085	if (activice) then
[633]	1086	pdt(1:ngrid,1:nlayer) =
	1087	& pdt(1:ngrid,1:nlayer) +
	1088	& zdtcloud(1:ngrid,1:nlayer)
[42]	1089	endif
[633]	1090
[42]	1091	! increment water vapour and ice atmospheric tracers tendencies
[706]	1092	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
	1093	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
	1094	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
	1095	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
	1096	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
	1097	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
	1098
	1099	! increment dust and ccn masses and numbers
[883]	1100	! We need to check that we have Nccn & Ndust > 0
	1101	! This is due to single precision rounding problems
[706]	1102	if (microphys) then
	1103	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
	1104	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
	1105	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
	1106	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
	1107	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
	1108	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
[883]	1109	where (pq(:,:,igcm_ccn_mass) +
	1110	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
	1111	pdq(:,:,igcm_ccn_mass) =
	1112	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
	1113	pdq(:,:,igcm_ccn_number) =
	1114	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1115	end where
	1116	where (pq(:,:,igcm_ccn_number) +
	1117	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
	1118	pdq(:,:,igcm_ccn_mass) =
	1119	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
	1120	pdq(:,:,igcm_ccn_number) =
	1121	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
	1122	end where
[706]	1123	endif
	1124
[883]	1125	if (scavenging) then
	1126	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
	1127	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
	1128	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
	1129	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
	1130	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
	1131	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
	1132	where (pq(:,:,igcm_dust_mass) +
	1133	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
	1134	pdq(:,:,igcm_dust_mass) =
	1135	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
	1136	pdq(:,:,igcm_dust_number) =
	1137	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1138	end where
	1139	where (pq(:,:,igcm_dust_number) +
	1140	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
	1141	pdq(:,:,igcm_dust_mass) =
	1142	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
	1143	pdq(:,:,igcm_dust_number) =
	1144	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
	1145	end where
	1146	endif ! of if scavenging
[633]	1147
	1148
[42]	1149	END IF ! of IF (water)
	1150
[556]	1151	c 7b. Aerosol particles
[42]	1152	c -------------------
	1153
	1154	c ----------
	1155	c Dust devil :
	1156	c ----------
	1157	IF(callddevil) then
[883]	1158	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
[42]	1159	& zdqdev,zdqsdev)
	1160
	1161	if (dustbin.ge.1) then
	1162	do iq=1,nq
	1163	DO l=1,nlayer
	1164	DO ig=1,ngrid
	1165	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	1166	ENDDO
	1167	ENDDO
	1168	enddo
	1169	do iq=1,nq
	1170	DO ig=1,ngrid
	1171	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
	1172	ENDDO
	1173	enddo
	1174	endif ! of if (dustbin.ge.1)
	1175
	1176	END IF ! of IF (callddevil)
	1177
	1178	c -------------
	1179	c Sedimentation : acts also on water ice
	1180	c -------------
	1181	IF (sedimentation) THEN
	1182	!call zerophys(ngridnlayernq, zdqsed)
	1183	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	1184	!call zerophys(ngrid*nq, zdqssed)
	1185	zdqssed(1:ngrid,1:nq)=0
	1186
	1187	call callsedim(ngrid,nlayer, ptimestep,
[1005]	1188	& zplev,zzlev, zzlay, pt, pdt, rdust, rice,
[358]	1189	& rsedcloud,rhocloud,
[411]	1190	& pq, pdq, zdqsed, zdqssed,nq,
	1191	& tau,tauscaling)
[833]	1192
[42]	1193	DO iq=1, nq
	1194	DO l=1,nlayer
	1195	DO ig=1,ngrid
	1196	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1197	ENDDO
	1198	ENDDO
	1199	ENDDO
	1200	DO iq=1, nq
	1201	DO ig=1,ngrid
	1202	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1203	ENDDO
	1204	ENDDO
	1205	END IF ! of IF (sedimentation)
[790]	1206
	1207	c Add lifted dust to tendancies after sedimentation in the LES
	1208	#ifdef MESOSCALE
	1209	#include "meso_inc/meso_inc_lift_les.F"
	1210	#endif
	1211
[556]	1212	c
	1213	c 7c. Chemical species
	1214	c ------------------
[42]	1215
[556]	1216	#ifndef MESOSCALE
	1217	c --------------
	1218	c photochemistry :
	1219	c --------------
	1220	IF (photochem .or. thermochem) then
[411]	1221
[556]	1222	! dust and ice surface area
[883]	1223	call surfacearea(ngrid, nlayer, ptimestep, zplay, zzlay,
[556]	1224	$ pt, pq, pdq, nq,
	1225	$ rdust, rice, tau, tauscaling,
	1226	$ surfdust, surfice)
	1227	! call photochemistry
[1036]	1228	call calchim(nq,
	1229	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
[556]	1230	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
	1231	$ zdqcloud,zdqscloud,tauref,co2ice,
	1232	$ pu,pdu,pv,pdv,surfdust,surfice)
	1233
	1234	! increment values of tracers:
	1235	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1236	! tracers is zero anyways
	1237	DO l=1,nlayer
	1238	DO ig=1,ngrid
	1239	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	1240	ENDDO
	1241	ENDDO
	1242	ENDDO ! of DO iq=1,nq
	1243
	1244	! add condensation tendency for H2O2
	1245	if (igcm_h2o2.ne.0) then
	1246	DO l=1,nlayer
	1247	DO ig=1,ngrid
	1248	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	1249	& +zdqcloud(ig,l,igcm_h2o2)
	1250	ENDDO
	1251	ENDDO
	1252	endif
	1253
	1254	! increment surface values of tracers:
	1255	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1256	! tracers is zero anyways
	1257	DO ig=1,ngrid
	1258	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	1259	ENDDO
	1260	ENDDO ! of DO iq=1,nq
	1261
	1262	! add condensation tendency for H2O2
	1263	if (igcm_h2o2.ne.0) then
	1264	DO ig=1,ngrid
	1265	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	1266	& +zdqscloud(ig,igcm_h2o2)
	1267	ENDDO
	1268	endif
	1269
	1270	END IF ! of IF (photochem.or.thermochem)
	1271	#endif
	1272
[42]	1273	c 7d. Updates
	1274	c ---------
	1275
	1276	DO iq=1, nq
	1277	DO ig=1,ngrid
	1278
	1279	c ---------------------------------
	1280	c Updating tracer budget on surface
	1281	c ---------------------------------
	1282	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1283
	1284	ENDDO ! (ig)
	1285	ENDDO ! (iq)
	1286
	1287	endif ! of if (tracer)
	1288
[226]	1289	#ifndef MESOSCALE
	1290	c-----------------------------------------------------------------------
	1291	c 8. THERMOSPHERE CALCULATION
	1292	c-----------------------------------------------------------------------
[42]	1293
[226]	1294	if (callthermos) then
[883]	1295	call thermosphere(zplev,zplay,dist_sol,
[226]	1296	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	1297	& pt,pq,pu,pv,pdt,pdq,
	1298	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	1299
	1300	DO l=1,nlayer
	1301	DO ig=1,ngrid
	1302	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	1303	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	1304	& +zdteuv(ig,l)
	1305	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	1306	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	1307	DO iq=1, nq
	1308	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	1309	ENDDO
	1310	ENDDO
	1311	ENDDO
	1312
	1313	endif ! of if (callthermos)
	1314	#endif
	1315
[42]	1316	c-----------------------------------------------------------------------
	1317	c 9. Surface and sub-surface soil temperature
	1318	c-----------------------------------------------------------------------
	1319	c
	1320	c
	1321	c 9.1 Increment Surface temperature:
	1322	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1323
	1324	DO ig=1,ngrid
	1325	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1326	ENDDO
	1327
	1328	c Prescribe a cold trap at south pole (except at high obliquity !!)
	1329	c Temperature at the surface is set there to be the temperature
	1330	c corresponding to equilibrium temperature between phases of CO2
	1331
[411]	1332
[42]	1333	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
[226]	1334	#ifndef MESOSCALE
	1335	if (caps.and.(obliquit.lt.27.)) then
	1336	! NB: Updated surface pressure, at grid point 'ngrid', is
[883]	1337	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
[226]	1338	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
[883]	1339	& (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
[226]	1340	endif
	1341	#endif
[42]	1342	c -------------------------------------------------------------
[283]	1343	c Change of surface albedo in case of ground frost
[42]	1344	c everywhere except on the north permanent cap and in regions
	1345	c covered by dry ice.
	1346	c ALWAYS PLACE these lines after newcondens !!!
	1347	c -------------------------------------------------------------
	1348	do ig=1,ngrid
	1349	if ((co2ice(ig).eq.0).and.
[283]	1350	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
	1351	albedo(ig,1) = albedo_h2o_ice
	1352	albedo(ig,2) = albedo_h2o_ice
	1353	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
	1354	c write(,) "physiq.F frost :"
	1355	c & ,lati(ig)180./pi, long(ig)180./pi
[42]	1356	endif
	1357	enddo ! of do ig=1,ngrid
	1358	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
	1359
[528]	1360	c
[42]	1361	c 9.2 Compute soil temperatures and subsurface heat flux:
	1362	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1363	IF (callsoil) THEN
[833]	1364	c Thermal inertia feedback
	1365	IF (tifeedback) THEN
	1366	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
	1367	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
	1368	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1369	ELSE
[42]	1370	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
[833]	1371	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1372	ENDIF
[42]	1373	ENDIF
[411]	1374
[42]	1375
[999]	1376
[42]	1377	c-----------------------------------------------------------------------
	1378	c 10. Write output files
	1379	c ----------------------
	1380
[674]	1381	c Save variables for eventual restart in MMM and LES
	1382	#ifdef MESOSCALE
	1383	#include "meso_inc/meso_inc_save_restart.F"
	1384	#endif
	1385
[42]	1386	c -------------------------------
	1387	c Dynamical fields incrementation
	1388	c -------------------------------
	1389	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1390	! temperature, zonal and meridional wind
	1391	DO l=1,nlayer
	1392	DO ig=1,ngrid
	1393	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1394	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1395	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1396	ENDDO
	1397	ENDDO
	1398
	1399	! tracers
	1400	DO iq=1, nq
	1401	DO l=1,nlayer
	1402	DO ig=1,ngrid
	1403	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1404	ENDDO
	1405	ENDDO
	1406	ENDDO
	1407
	1408	! Density
	1409	DO l=1,nlayer
	1410	DO ig=1,ngrid
	1411	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1412	ENDDO
	1413	ENDDO
	1414
[269]	1415	! Potential Temperature
	1416
	1417	DO ig=1,ngridmx
	1418	DO l=1,nlayermx
[528]	1419	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
[269]	1420	ENDDO
	1421	ENDDO
	1422
	1423
[226]	1424	c Compute surface stress : (NB: z0 is a common in surfdat.h)
[42]	1425	c DO ig=1,ngrid
[226]	1426	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
[42]	1427	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1428	c ENDDO
	1429
	1430	c Sum of fluxes in solar spectral bands (for output only)
	1431	DO ig=1,ngrid
[226]	1432	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1433	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
[42]	1434	ENDDO
	1435	c ***** TEST ****************************************************
[835]	1436	ztim1 = 999
	1437	DO l=1,nlayer
	1438	DO ig=1,ngrid
	1439	if (pt(ig,l).lt.ztim1) then
	1440	ztim1 = pt(ig,l)
	1441	igmin = ig
	1442	lmin = l
	1443	end if
	1444	ENDDO
	1445	ENDDO
	1446	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
	1447	write(,) 'PHYSIQ: stability WARNING :'
	1448	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1449	& 'ig l =', igmin, lmin
	1450	end if
[42]	1451	c *******************************************************************
	1452
	1453	c ---------------------
	1454	c Outputs to the screen
	1455	c ---------------------
	1456
	1457	IF (lwrite) THEN
	1458	PRINT*,'Global diagnostics for the physics'
	1459	PRINT,'Variables and their increments x and dx/dt dt'
	1460	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1461	WRITE(*,'(2f10.5,2f15.5)')
	1462	s tsurf(igout),zdtsurf(igout)*ptimestep,
[883]	1463	s zplev(igout,1),pdpsrf(igout)*ptimestep
[42]	1464	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1465	WRITE(*,'(i4,6f10.5)') (l,
	1466	s pu(igout,l),pdu(igout,l)*ptimestep,
	1467	s pv(igout,l),pdv(igout,l)*ptimestep,
	1468	s pt(igout,l),pdt(igout,l)*ptimestep,
	1469	s l=1,nlayer)
	1470	ENDIF ! of IF (lwrite)
	1471
[566]	1472	c ----------------------------------------------------------
	1473	c ----------------------------------------------------------
	1474	c INTERPOLATIONS IN THE SURFACE-LAYER
	1475	c ----------------------------------------------------------
	1476	c ----------------------------------------------------------
	1477
[657]	1478	n_out=0 ! number of elements in the z_out array.
[636]	1479	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
	1480	! to 5
	1481	IF (n_out .ne. 0) THEN
[566]	1482
[636]	1483	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
	1484	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
	1485	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
[566]	1486
[636]	1487	z_out(:)=[3.,2.,1.,0.5,0.1]
	1488	u_out(:,:)=0.
	1489	T_out(:,:)=0.
[566]	1490
[636]	1491	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
[566]	1492	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
[657]	1493	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
[566]	1494
	1495	#ifndef MESOSCALE
	1496	IF (ngrid .eq. 1) THEN
	1497	dimout=0
	1498	ELSE
	1499	dimout=2
	1500	ENDIF
	1501	DO n=1,n_out
[636]	1502	write(zstring, '(F8.6)') z_out(n)
	1503	call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
	1504	& 'potential temperature at z_out','K',dimout,T_out(:,n))
[566]	1505	call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
	1506	& 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
[636]	1507	ENDDO
[566]	1508	call WRITEDIAGFI(ngrid,'u_star',
	1509	& 'friction velocity','m/s',dimout,ustar)
	1510	call WRITEDIAGFI(ngrid,'teta_star',
	1511	& 'friction potential temperature','K',dimout,tstar)
	1512	call WRITEDIAGFI(ngrid,'L',
	1513	& 'Monin Obukhov length','m',dimout,L_mo)
	1514	call WRITEDIAGFI(ngrid,'vvv',
	1515	& 'Vertical velocity variance at zout','m',dimout,vvv)
	1516	call WRITEDIAGFI(ngrid,'vhf',
	1517	& 'Vertical heat flux at zout','m',dimout,vhf)
[636]	1518	#else
	1519	T_out1(:)=T_out(:,1)
	1520	u_out1(:)=u_out(:,1)
[566]	1521	#endif
	1522
	1523	ENDIF
	1524
	1525	c ----------------------------------------------------------
	1526	c ----------------------------------------------------------
	1527	c END OF SURFACE LAYER INTERPOLATIONS
	1528	c ----------------------------------------------------------
	1529	c ----------------------------------------------------------
	1530
[42]	1531	IF (ngrid.NE.1) THEN
	1532
[226]	1533	#ifndef MESOSCALE
[42]	1534	c -------------------------------------------------------------------
	1535	c Writing NetCDF file "RESTARTFI" at the end of the run
	1536	c -------------------------------------------------------------------
	1537	c Note: 'restartfi' is stored just before dynamics are stored
	1538	c in 'restart'. Between now and the writting of 'restart',
	1539	c there will have been the itau=itau+1 instruction and
	1540	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1541	c thus we store for time=time+dtvr
	1542
[999]	1543	IF( ((ecritstart.GT.0) .and.
	1544	. (MOD(icount*iphysiq,ecritstart).EQ.0))
	1545	. .or. lastcall ) THEN
	1546
	1547	ztime_fin = pday + ptime + ptimestep/(float(iphysiq)*daysec)
	1548	. - day_ini - time_phys
	1549	print*, pday,ptime,day_ini, time_phys
	1550	write(*,'(A,I7,A,F12.5)')
	1551	. 'PHYSIQ: Ecriture du fichier restartfi ; icount=',
	1552	. icount,' date=',ztime_fin
	1553
	1554
	1555	call physdem1("restartfi.nc",nsoilmx,nq,
	1556	. ptimestep,ztime_fin,
	1557	. tsurf,tsoil,co2ice,emis,q2,qsurf)
	1558
[226]	1559	ENDIF
	1560	#endif
[42]	1561
	1562	c -------------------------------------------------------------------
	1563	c Calculation of diagnostic variables written in both stats and
	1564	c diagfi files
	1565	c -------------------------------------------------------------------
	1566
	1567	if (tracer) then
[719]	1568
	1569	if(doubleq) then
	1570	do ig=1,ngrid
	1571	dqdustsurf(ig) =
[756]	1572	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
[719]	1573	dndustsurf(ig) =
[756]	1574	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
	1575	ndust(ig,:) =
	1576	& pq(ig,:,igcm_dust_number)*tauscaling(ig)
	1577	qdust(ig,:) =
	1578	& pq(ig,:,igcm_dust_mass)*tauscaling(ig)
[719]	1579	enddo
	1580	if (scavenging) then
	1581	do ig=1,ngrid
	1582	dqdustsurf(ig) = dqdustsurf(ig) +
[756]	1583	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
[719]	1584	dndustsurf(ig) = dndustsurf(ig) +
[756]	1585	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
	1586	nccn(ig,:) =
	1587	& pq(ig,:,igcm_ccn_number)*tauscaling(ig)
	1588	qccn(ig,:) =
	1589	& pq(ig,:,igcm_ccn_mass)*tauscaling(ig)
[719]	1590	enddo
	1591	endif
	1592	endif
	1593
[42]	1594	if (water) then
[286]	1595	mtot(:)=0
	1596	icetot(:)=0
	1597	rave(:)=0
	1598	tauTES(:)=0
[42]	1599	do ig=1,ngrid
	1600	do l=1,nlayermx
	1601	mtot(ig) = mtot(ig) +
	1602	& zq(ig,l,igcm_h2o_vap) *
[883]	1603	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	1604	icetot(ig) = icetot(ig) +
	1605	& zq(ig,l,igcm_h2o_ice) *
[883]	1606	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	1607	c Computing abs optical depth at 825 cm-1 in each
	1608	c layer to simulate NEW TES retrieval
	1609	Qabsice = min(
	1610	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1611	& )
	1612	opTES(ig,l)= 0.75 * Qabsice *
	1613	& zq(ig,l,igcm_h2o_ice) *
[883]	1614	& (zplev(ig,l) - zplev(ig,l+1)) / g
[42]	1615	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1616	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1617	enddo
[669]	1618	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
	1619	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
[42]	1620	enddo
[883]	1621	call watersat(ngridmx*nlayermx,zt,zplay,zqsat)
[669]	1622	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
[42]	1623
[669]	1624	if (scavenging) then
	1625	Nccntot(:)= 0
	1626	Mccntot(:)= 0
	1627	rave(:)=0
	1628	do ig=1,ngrid
	1629	do l=1,nlayermx
	1630	Nccntot(ig) = Nccntot(ig) +
	1631	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
[883]	1632	& *(zplev(ig,l) - zplev(ig,l+1)) / g
[669]	1633	Mccntot(ig) = Mccntot(ig) +
	1634	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
[883]	1635	& *(zplev(ig,l) - zplev(ig,l+1)) / g
[669]	1636	cccc Column integrated effective ice radius
	1637	cccc is weighted by total ice surface area (BETTER than total ice mass)
	1638	rave(ig) = rave(ig) +
	1639	& tauscaling(ig) *
	1640	& zq(ig,l,igcm_ccn_number) *
[883]	1641	& (zplev(ig,l) - zplev(ig,l+1)) / g *
[669]	1642	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
	1643	enddo
	1644	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
	1645	& /max(pi*rave(ig),1.e-30) ! surface weight
	1646	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1647	enddo
[833]	1648	else ! of if (scavenging)
	1649	rave(:)=0
	1650	do ig=1,ngrid
	1651	do l=1,nlayermx
	1652	rave(ig) = rave(ig) +
	1653	& zq(ig,l,igcm_h2o_ice) *
[883]	1654	& (zplev(ig,l) - zplev(ig,l+1)) / g *
[833]	1655	& rice(ig,l) * (1.+nuice_ref)
	1656	enddo
	1657	rave(ig) = max(rave(ig) /
	1658	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
	1659	enddo
[669]	1660	endif ! of if (scavenging)
	1661
[42]	1662	endif ! of if (water)
	1663	endif ! of if (tracer)
	1664
	1665	c -----------------------------------------------------------------
	1666	c WSTATS: Saving statistics
	1667	c -----------------------------------------------------------------
	1668	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1669	c which can later be used to make the statistic files of the run:
	1670	c "stats") only possible in 3D runs !
	1671
[695]	1672	IF (callstats) THEN
[42]	1673
[695]	1674	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1675	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1676	call wstats(ngrid,"co2ice","CO2 ice cover",
[575]	1677	& "kg.m-2",2,co2ice)
[758]	1678	call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
	1679	& 2,tauref)
[695]	1680	call wstats(ngrid,"fluxsurf_lw",
[575]	1681	& "Thermal IR radiative flux to surface","W.m-2",2,
	1682	& fluxsurf_lw)
[695]	1683	call wstats(ngrid,"fluxsurf_sw",
[575]	1684	& "Solar radiative flux to surface","W.m-2",2,
	1685	& fluxsurf_sw_tot)
[695]	1686	call wstats(ngrid,"fluxtop_lw",
[575]	1687	& "Thermal IR radiative flux to space","W.m-2",2,
	1688	& fluxtop_lw)
[695]	1689	call wstats(ngrid,"fluxtop_sw",
[575]	1690	& "Solar radiative flux to space","W.m-2",2,
	1691	& fluxtop_sw_tot)
[695]	1692	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1693	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1694	call wstats(ngrid,"v","Meridional (North-South) wind",
[575]	1695	& "m.s-1",3,zv)
[705]	1696	call wstats(ngrid,"w","Vertical (down-up) wind",
	1697	& "m.s-1",3,pw)
[695]	1698	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
[883]	1699	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
[758]	1700	call wstats(ngrid,"q2",
	1701	& "Boundary layer eddy kinetic energy",
	1702	& "m2.s-2",3,q2)
	1703	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1704	& emis)
[226]	1705	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1706	c & 2,zstress)
	1707	c call wstats(ngrid,"sw_htrt","sw heat.rate",
	1708	c & "W.m-2",3,zdtsw)
	1709	c call wstats(ngrid,"lw_htrt","lw heat.rate",
	1710	c & "W.m-2",3,zdtlw)
[42]	1711
[758]	1712	if (calltherm) then
	1713	call wstats(ngrid,"zmax_th","Height of thermals",
	1714	& "m",2,zmax_th)
	1715	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
	1716	& "K.m/s",2,hfmax_th)
	1717	call wstats(ngrid,"wstar",
	1718	& "Max vertical velocity in thermals",
	1719	& "m/s",2,wstar)
	1720	endif
	1721
[226]	1722	if (tracer) then
	1723	if (water) then
[520]	1724	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
[719]	1725	& *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_vap)
	1726	call wstats(ngrid,"vmr_h2ovap",
[520]	1727	& "H2O vapor volume mixing ratio","mol/mol",
	1728	& 3,vmr)
	1729	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
[719]	1730	& *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_ice)
[520]	1731	call wstats(ngrid,"vmr_h2oice",
	1732	& "H2O ice volume mixing ratio","mol/mol",
	1733	& 3,vmr)
[669]	1734	vmr=zqsat(1:ngridmx,1:nlayermx)
[719]	1735	& *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_vap)
[669]	1736	call wstats(ngrid,"vmr_h2osat",
	1737	& "saturation volume mixing ratio","mol/mol",
	1738	& 3,vmr)
[226]	1739	call wstats(ngrid,"h2o_ice_s",
	1740	& "surface h2o_ice","kg/m2",
	1741	& 2,qsurf(1,igcm_h2o_ice))
[758]	1742	call wstats(ngrid,'albedo',
	1743	& 'albedo',
[1013]	1744	& '',2,albedo(1,1))
[226]	1745	call wstats(ngrid,"mtot",
	1746	& "total mass of water vapor","kg/m2",
	1747	& 2,mtot)
	1748	call wstats(ngrid,"icetot",
	1749	& "total mass of water ice","kg/m2",
	1750	& 2,icetot)
[520]	1751	call wstats(ngrid,"reffice",
	1752	& "Mean reff","m",
	1753	& 2,rave)
[669]	1754	call wstats(ngrid,"Nccntot",
[520]	1755	& "condensation nuclei","Nbr/m2",
[669]	1756	& 2,Nccntot)
	1757	call wstats(ngrid,"Mccntot",
	1758	& "condensation nuclei mass","kg/m2",
	1759	& 2,Mccntot)
[520]	1760	call wstats(ngrid,"rice",
	1761	& "Ice particle size","m",
	1762	& 3,rice)
[226]	1763	if (.not.activice) then
	1764	call wstats(ngrid,"tauTESap",
	1765	& "tau abs 825 cm-1","",
	1766	& 2,tauTES)
[520]	1767	else
	1768	call wstats(ngridmx,'tauTES',
	1769	& 'tau abs 825 cm-1',
	1770	& '',2,taucloudtes)
[226]	1771	endif
	1772
	1773	endif ! of if (water)
[719]	1774
	1775
	1776	if (dustbin.ne.0) then
	1777
	1778	call wstats(ngridmx,'tau','taudust','SI',2,tau(1,1))
	1779
	1780	if (doubleq) then
	1781	c call wstats(ngridmx,'qsurf','qsurf',
	1782	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
	1783	c call wstats(ngridmx,'Nsurf','N particles',
	1784	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
	1785	c call wstats(ngridmx,'dqsdev','ddevil lift',
	1786	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
	1787	c call wstats(ngridmx,'dqssed','sedimentation',
	1788	c & 'kg.m-2.s-1',2,zdqssed(1,1))
	1789	c call wstats(ngridmx,'dqsdif','diffusion',
	1790	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
	1791	call wstats(ngridmx,'dqsdust',
	1792	& 'deposited surface dust mass',
	1793	& 'kg.m-2.s-1',2,dqdustsurf)
	1794	call wstats(ngridmx,'dqndust',
	1795	& 'deposited surface dust number',
	1796	& 'number.m-2.s-1',2,dndustsurf)
	1797	call wstats(ngridmx,'reffdust','reffdust',
	1798	& 'm',3,rdust*ref_r0)
	1799	call wstats(ngridmx,'dustq','Dust mass mr',
[756]	1800	& 'kg/kg',3,qdust)
[719]	1801	call wstats(ngridmx,'dustN','Dust number',
[756]	1802	& 'part/kg',3,ndust)
[719]	1803	else
	1804	do iq=1,dustbin
	1805	write(str2(1:2),'(i2.2)') iq
	1806	call wstats(ngridmx,'q'//str2,'mix. ratio',
	1807	& 'kg/kg',3,zq(1,1,iq))
	1808	call wstats(ngridmx,'qsurf'//str2,'qsurf',
	1809	& 'kg.m-2',2,qsurf(1,iq))
	1810	end do
	1811	endif ! (doubleq)
[226]	1812
[719]	1813	if (scavenging) then
	1814	call wstats(ngridmx,'ccnq','CCN mass mr',
[756]	1815	& 'kg/kg',3,qccn)
[719]	1816	call wstats(ngridmx,'ccnN','CCN number',
[756]	1817	& 'part/kg',3,nccn)
[719]	1818	endif ! (scavenging)
	1819
	1820	endif ! (dustbin.ne.0)
	1821
	1822
[226]	1823	if (thermochem.or.photochem) then
	1824	do iq=1,nq
[576]	1825	if (noms(iq) .ne. "dust_mass" .and.
	1826	$ noms(iq) .ne. "dust_number" .and.
	1827	$ noms(iq) .ne. "ccn_mass" .and.
[719]	1828	$ noms(iq) .ne. "ccn_number" .and.
	1829	$ noms(iq) .ne. "h2o_vap" .and.
	1830	$ noms(iq) .ne. "h2o_ice") then
[705]	1831	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	1832	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
	1833	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
	1834	& rho(1:ngrid,1:nlayer)n_avog/
	1835	& (1000*mmol(iq))
[576]	1836	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1837	$ "Volume mixing ratio","mol/mol",3,vmr)
[705]	1838	! call wstats(ngrid,"rho_"//trim(noms(iq)),
	1839	! $ "Number density","cm-3",3,rhopart)
	1840	! call writediagfi(ngrid,"rho_"//trim(noms(iq)),
	1841	! $ "Number density","cm-3",3,rhopart)
[576]	1842	if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
	1843	$ (noms(iq).eq."o3")) then
	1844	call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
	1845	$ "Volume mixing ratio","mol/mol",3,vmr)
	1846	end if
	1847	do ig = 1,ngrid
	1848	colden(ig,iq) = 0.
	1849	end do
	1850	do l=1,nlayer
	1851	do ig=1,ngrid
	1852	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
[883]	1853	$ *(zplev(ig,l)-zplev(ig,l+1))
[576]	1854	$ 6.022e22/(mmol(iq)g)
	1855	end do
	1856	end do
	1857	call wstats(ngrid,"c_"//trim(noms(iq)),
	1858	$ "column","mol cm-2",2,colden(1,iq))
	1859	call writediagfi(ngrid,"c_"//trim(noms(iq)),
	1860	$ "column","mol cm-2",2,colden(1,iq))
	1861	end if ! of if (noms(iq) .ne. "dust_mass" ...)
	1862	end do ! of do iq=1,nq
[334]	1863	end if ! of if (thermochem.or.photochem)
[226]	1864
[334]	1865	end if ! of if (tracer)
[226]	1866
	1867	IF(lastcall) THEN
	1868	write (,) "Writing stats..."
	1869	call mkstats(ierr)
	1870	ENDIF
	1871
[42]	1872	ENDIF !if callstats
	1873
	1874	c (Store EOF for Mars Climate database software)
	1875	IF (calleofdump) THEN
	1876	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1877	ENDIF
	1878
[234]	1879
[226]	1880	#ifdef MESOSCALE
[234]	1881	!!!
	1882	!!! OUTPUT FIELDS
	1883	!!!
[81]	1884	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
[668]	1885	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
	1886	TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref)
[667]	1887	IF (tracer) THEN
[82]	1888	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
[608]	1889	icetot(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
[308]	1890	!! JF
	1891	IF (igcm_dust_mass .ne. 0) THEN
[708]	1892	qsurfdust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
[308]	1893	ENDIF
[81]	1894	IF (igcm_h2o_ice .ne. 0) THEN
	1895	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
	1896	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
[719]	1897	. *mmean(1:ngridmx,1:nlayermx) / mmol(igcm_h2o_ice)
[81]	1898	ENDIF
[308]	1899	!! Dust quantity integration along the vertical axe
	1900	dustot(:)=0
[698]	1901	IF (igcm_dust_mass .ne. 0) THEN
[308]	1902	do ig=1,ngrid
	1903	do l=1,nlayermx
	1904	dustot(ig) = dustot(ig) +
	1905	& zq(ig,l,igcm_dust_mass)
[883]	1906	& * (zplev(ig,l) - zplev(ig,l+1)) / g
[308]	1907	enddo
	1908	enddo
[667]	1909	ENDIF
[698]	1910	ENDIF
[550]	1911	!! TAU water ice as seen by TES
	1912	if (activice) tauTES = taucloudtes
[226]	1913	c AUTOMATICALLY GENERATED FROM REGISTRY
[42]	1914	#include "fill_save.inc"
[226]	1915	#else
[528]	1916	#ifndef MESOINI
[42]	1917
[226]	1918	c ==========================================================
	1919	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1920	c any variable for diagnostic (output with period
	1921	c "ecritphy", set in "run.def")
	1922	c ==========================================================
	1923	c WRITEDIAGFI can ALSO be called from any other subroutines
	1924	c for any variables !!
[835]	1925	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1926	c & emis)
	1927	c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
	1928	c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
[226]	1929	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1930	& tsurf)
	1931	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
[769]	1932	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
	1933	& ,"kg.m-2",2,co2ice)
[284]	1934
[719]	1935	call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1936	& "K",2,zt(1,7))
	1937	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1938	& fluxsurf_lw)
	1939	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1940	& fluxsurf_sw_tot)
	1941	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1942	& fluxtop_lw)
	1943	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1944	& fluxtop_sw_tot)
	1945	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1946	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1947	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1948	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1949	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
[226]	1950	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
[520]	1951	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
[758]	1952	call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
[226]	1953	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1954	c & zstress)
[835]	1955	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	1956	c & 'w.m-2',3,zdtsw)
	1957	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	1958	c & 'w.m-2',3,zdtlw)
[520]	1959	if (.not.activice) then
	1960	CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1961	& 'tau abs 825 cm-1',
	1962	& '',2,tauTES)
	1963	else
	1964	CALL WRITEDIAGFI(ngridmx,'tauTES',
	1965	& 'tau abs 825 cm-1',
	1966	& '',2,taucloudtes)
	1967	endif
[42]	1968
[528]	1969	#else
	1970	!!! this is to ensure correct initialisation of mesoscale model
	1971	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1972	& tsurf)
	1973	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1974	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1975	& co2ice)
	1976	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1977	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1978	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
[299]	1979	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1980	& emis)
	1981	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
	1982	& "K",3,tsoil)
	1983	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
	1984	& "K",3,inertiedat)
	1985	#endif
[42]	1986
[299]	1987
[226]	1988	c ----------------------------------------------------------
	1989	c Outputs of the CO2 cycle
	1990	c ----------------------------------------------------------
[42]	1991
[226]	1992	if (tracer.and.(igcm_co2.ne.0)) then
	1993	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	1994	! & "kg/kg",2,zq(1,1,igcm_co2))
[719]	1995	call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	1996	& "kg/kg",3,zq(1,1,igcm_co2))
[226]	1997
	1998	! Compute co2 column
[286]	1999	co2col(:)=0
[226]	2000	do l=1,nlayermx
	2001	do ig=1,ngrid
	2002	co2col(ig)=co2col(ig)+
[883]	2003	& zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
[226]	2004	enddo
	2005	enddo
	2006	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	2007	& co2col)
	2008	endif ! of if (tracer.and.(igcm_co2.ne.0))
[42]	2009
[226]	2010	c ----------------------------------------------------------
	2011	c Outputs of the water cycle
	2012	c ----------------------------------------------------------
	2013	if (tracer) then
	2014	if (water) then
	2015
[299]	2016	#ifdef MESOINI
	2017	!!!! waterice = q01, voir readmeteo.F90
	2018	call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
	2019	& 'kg/kg',3,
	2020	& zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
	2021	!!!! watervapor = q02, voir readmeteo.F90
	2022	call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap),
	2023	& 'kg/kg',3,
	2024	& zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
	2025	!!!! surface waterice qsurf02 (voir readmeteo)
	2026	call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer',
	2027	& 'kg.m-2',2,
	2028	& qsurf(1:ngridmx,igcm_h2o_ice))
	2029	#endif
	2030
[719]	2031	CALL WRITEDIAGFI(ngridmx,'mtot',
[226]	2032	& 'total mass of water vapor',
	2033	& 'kg/m2',2,mtot)
[719]	2034	CALL WRITEDIAGFI(ngridmx,'icetot',
[226]	2035	& 'total mass of water ice',
	2036	& 'kg/m2',2,icetot)
[719]	2037	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	2038	& *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_ice)
	2039	call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
	2040	& 'mol/mol',3,vmr)
	2041	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	2042	& *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_vap)
	2043	call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
	2044	& 'mol/mol',3,vmr)
	2045	CALL WRITEDIAGFI(ngridmx,'reffice',
[669]	2046	& 'Mean reff',
	2047	& 'm',2,rave)
[719]	2048	CALL WRITEDIAGFI(ngrid,"Nccntot",
[669]	2049	& "condensation nuclei","Nbr/m2",
	2050	& 2,Nccntot)
[719]	2051	CALL WRITEDIAGFI(ngrid,"Mccntot",
	2052	& "mass condensation nuclei","kg/m2",
	2053	& 2,Mccntot)
	2054	call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
	2055	& 'm',3,rice)
	2056	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
[226]	2057	& 'surface h2o_ice',
	2058	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
[833]	2059	CALL WRITEDIAGFI(ngridmx,'albedo',
	2060	& 'albedo',
[1013]	2061	& '',2,albedo(1,1))
[833]	2062	if (tifeedback) then
	2063	call WRITEDIAGSOIL(ngridmx,"soiltemp",
	2064	& "Soil temperature","K",
	2065	& 3,tsoil)
	2066	call WRITEDIAGSOIL(ngridmx,'soilti',
	2067	& 'Soil Thermal Inertia',
	2068	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
	2069	endif
[226]	2070	endif !(water)
	2071
	2072
	2073	if (water.and..not.photochem) then
	2074	iq=nq
	2075	c write(str2(1:2),'(i2.2)') iq
	2076	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
	2077	c & 'kg.m-2',2,zdqscloud(1,iq))
	2078	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
	2079	c & 'kg/kg',3,zdqchim(1,1,iq))
	2080	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
	2081	c & 'kg/kg',3,zdqdif(1,1,iq))
	2082	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
	2083	c & 'kg/kg',3,zdqadj(1,1,iq))
	2084	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
	2085	c & 'kg/kg',3,zdqc(1,1,iq))
	2086	endif !(water.and..not.photochem)
	2087	endif
	2088
	2089	c ----------------------------------------------------------
	2090	c Outputs of the dust cycle
	2091	c ----------------------------------------------------------
	2092
[758]	2093	call WRITEDIAGFI(ngridmx,'tauref',
	2094	& 'Dust ref opt depth','NU',2,tauref)
[226]	2095
	2096	if (tracer.and.(dustbin.ne.0)) then
[719]	2097	call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
[226]	2098	if (doubleq) then
[411]	2099	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
	2100	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
	2101	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
	2102	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
[226]	2103	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
	2104	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
[411]	2105	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
	2106	c & 'kg.m-2.s-1',2,zdqssed(1,1))
	2107	c call WRITEDIAGFI(ngridmx,'dqsdif','diffusion',
	2108	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
[719]	2109	call WRITEDIAGFI(ngridmx,'dqsdust',
	2110	& 'deposited surface dust mass',
	2111	& 'kg.m-2.s-1',2,dqdustsurf)
	2112	call WRITEDIAGFI(ngridmx,'dqndust',
	2113	& 'deposited surface dust number',
	2114	& 'number.m-2.s-1',2,dndustsurf)
	2115	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
	2116	& 'm',3,rdust*ref_r0)
	2117	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
[756]	2118	& 'kg/kg',3,qdust)
[719]	2119	call WRITEDIAGFI(ngridmx,'dustN','Dust number',
[756]	2120	& 'part/kg',3,ndust)
[299]	2121	#ifdef MESOINI
[1032]	2122	! !!! to initialize mesoscale we need scaled variables
	2123	! !!! because this must correspond to starting point for tracers
	2124	! call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	2125	! & 'kg/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_mass))
	2126	! call WRITEDIAGFI(ngridmx,'dustN','Dust number',
	2127	! & 'part/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_number))
	2128	! call WRITEDIAGFI(ngridmx,'ccn','Nuclei mass mr',
	2129	! & 'kg/kg',3,pq(1:ngridmx,1:nlayermx,igcm_ccn_mass))
	2130	! call WRITEDIAGFI(ngridmx,'ccnN','Nuclei number',
	2131	! & 'part/kg',3,pq(1:ngridmx,1:nlayermx,igcm_ccn_number))
	2132	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	2133	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
	2134	call WRITEDIAGFI(ngridmx,'dustN','Dust number',
	2135	& 'part/kg',3,pq(1,1,igcm_dust_number))
	2136	call WRITEDIAGFI(ngridmx,'ccn','Nuclei mass mr',
	2137	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
	2138	call WRITEDIAGFI(ngridmx,'ccnN','Nuclei number',
	2139	& 'part/kg',3,pq(1,1,igcm_ccn_number))
[299]	2140	#endif
[226]	2141	else
	2142	do iq=1,dustbin
	2143	write(str2(1:2),'(i2.2)') iq
	2144	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
	2145	& 'kg/kg',3,zq(1,1,iq))
	2146	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
	2147	& 'kg.m-2',2,qsurf(1,iq))
	2148	end do
	2149	endif ! (doubleq)
[358]	2150
	2151	if (scavenging) then
[719]	2152	call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
[756]	2153	& 'kg/kg',3,qccn)
[719]	2154	call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
[756]	2155	& 'part/kg',3,nccn)
[358]	2156	endif ! (scavenging)
	2157
[226]	2158	c if (submicron) then
	2159	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
	2160	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
	2161	c endif ! (submicron)
	2162	end if ! (tracer.and.(dustbin.ne.0))
	2163
[705]	2164
[226]	2165	c ----------------------------------------------------------
[705]	2166	c Thermospheric outputs
[267]	2167	c ----------------------------------------------------------
[705]	2168
	2169	if(callthermos) then
	2170
	2171	call WRITEDIAGFI(ngridmx,"q15um","15 um cooling","K/s",
	2172	$ 3,zdtnlte)
	2173	call WRITEDIAGFI(ngridmx,"quv","UV heating","K/s",
	2174	$ 3,zdteuv)
	2175	call WRITEDIAGFI(ngridmx,"cond","Thermal conduction","K/s",
	2176	$ 3,zdtconduc)
	2177	call WRITEDIAGFI(ngridmx,"qnir","NIR heating","K/s",
	2178	$ 3,zdtnirco2)
	2179
	2180	endif !(callthermos)
	2181
	2182	c ----------------------------------------------------------
	2183	c ----------------------------------------------------------
[267]	2184	c PBL OUTPUS
	2185	c ----------------------------------------------------------
	2186	c ----------------------------------------------------------
	2187
	2188	c ----------------------------------------------------------
[226]	2189	c Outputs of thermals
	2190	c ----------------------------------------------------------
	2191	if (calltherm) then
	2192
	2193	! call WRITEDIAGFI(ngrid,'dtke',
	2194	! & 'tendance tke thermiques','m2/s2',
	2195	! & 3,dtke_th)
	2196	! call WRITEDIAGFI(ngrid,'d_u_ajs',
	2197	! & 'tendance u thermiques','m/s',
	2198	! & 3,pdu_th*ptimestep)
	2199	! call WRITEDIAGFI(ngrid,'d_v_ajs',
	2200	! & 'tendance v thermiques','m/s',
	2201	! & 3,pdv_th*ptimestep)
	2202	! if (tracer) then
	2203	! if (nq .eq. 2) then
	2204	! call WRITEDIAGFI(ngrid,'deltaq_th',
	2205	! & 'delta q thermiques','kg/kg',
	2206	! & 3,ptimestep*pdq_th(:,:,2))
	2207	! endif
	2208	! endif
	2209
[523]	2210	call WRITEDIAGFI(ngridmx,'zmax_th',
	2211	& 'hauteur du thermique','m',
	2212	& 2,zmax_th)
[226]	2213	call WRITEDIAGFI(ngridmx,'hfmax_th',
	2214	& 'maximum TH heat flux','K.m/s',
	2215	& 2,hfmax_th)
[528]	2216	call WRITEDIAGFI(ngridmx,'wstar',
[226]	2217	& 'maximum TH vertical velocity','m/s',
[528]	2218	& 2,wstar)
[226]	2219
	2220	endif
	2221
	2222	c ----------------------------------------------------------
[267]	2223	c ----------------------------------------------------------
	2224	c END OF PBL OUTPUS
	2225	c ----------------------------------------------------------
	2226	c ----------------------------------------------------------
	2227
	2228
	2229	c ----------------------------------------------------------
[226]	2230	c Output in netcdf file "diagsoil.nc" for subterranean
	2231	c variables (output every "ecritphy", as for writediagfi)
	2232	c ----------------------------------------------------------
	2233
	2234	! Write soil temperature
	2235	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	2236	! & 3,tsoil)
	2237	! Write surface temperature
	2238	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	2239	! & 2,tsurf)
	2240
	2241	c ==========================================================
	2242	c END OF WRITEDIAGFI
	2243	c ==========================================================
	2244	#endif
	2245
[42]	2246	ELSE ! if(ngrid.eq.1)
	2247
[627]	2248	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
	2249	& zls*180./pi,odpref,tauref
[42]	2250	c ----------------------------------------------------------------------
	2251	c Output in grads file "g1d" (ONLY when using testphys1d)
	2252	c (output at every X physical timestep)
	2253	c ----------------------------------------------------------------------
	2254	c
	2255	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	2256	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	2257	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	2258
	2259	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	2260	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	2261	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	2262	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	2263
[226]	2264	! THERMALS STUFF 1D
	2265	if(calltherm) then
	2266
	2267	call WRITEDIAGFI(ngridmx,'lmax_th',
	2268	& 'hauteur du thermique','point',
	2269	& 0,lmax_th_out)
[528]	2270	call WRITEDIAGFI(ngridmx,'zmax_th',
	2271	& 'hauteur du thermique','m',
	2272	& 0,zmax_th)
[226]	2273	call WRITEDIAGFI(ngridmx,'hfmax_th',
	2274	& 'maximum TH heat flux','K.m/s',
	2275	& 0,hfmax_th)
[528]	2276	call WRITEDIAGFI(ngridmx,'wstar',
[226]	2277	& 'maximum TH vertical velocity','m/s',
[528]	2278	& 0,wstar)
[226]	2279
	2280	co2col(:)=0.
	2281	if (tracer) then
	2282	do l=1,nlayermx
	2283	do ig=1,ngrid
	2284	co2col(ig)=co2col(ig)+
[883]	2285	& zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
[226]	2286	enddo
	2287	enddo
	2288
	2289	end if
	2290	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
	2291	& ,'kg/m-2',0,co2col)
[722]	2292	endif ! of if (calltherm)
	2293
[226]	2294	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
	2295	& ,'m/s',1,pw)
	2296	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
	2297	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
	2298	& tsurf)
[277]	2299	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
	2300	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
[226]	2301
	2302	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
	2303	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
[358]	2304	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
[544]	2305	! call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
	2306	! & "K.s-1",1,dtrad/zpopsk)
	2307	! call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	2308	! & 'w.m-2',1,zdtsw/zpopsk)
	2309	! call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	2310	! & 'w.m-2',1,zdtlw/zpopsk)
[769]	2311	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
	2312	& ,"kg.m-2",0,co2ice)
[544]	2313
[226]	2314	! or output in diagfi.nc (for testphys1d)
	2315	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
	2316	call WRITEDIAGFI(ngridmx,'temp','Temperature',
	2317	& 'K',1,zt)
[358]	2318
[226]	2319	if(tracer) then
	2320	c CALL writeg1d(ngrid,1,tau,'tau','SI')
	2321	do iq=1,nq
	2322	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	2323	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
	2324	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	2325	end do
[358]	2326	if (doubleq) then
	2327	call WRITEDIAGFI(ngridmx,'rdust','rdust',
	2328	& 'm',1,rdust)
	2329	endif
[833]	2330	if (water.AND.tifeedback) then
	2331	call WRITEDIAGFI(ngridmx,"soiltemp",
	2332	& "Soil temperature","K",
	2333	& 1,tsoil)
	2334	call WRITEDIAGFI(ngridmx,'soilti',
	2335	& 'Soil Thermal Inertia',
	2336	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
	2337	endif
[226]	2338	end if
[358]	2339
[520]	2340	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
[358]	2341	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[520]	2342	IF (water) THEN
[722]	2343
	2344	if (.not.activice) then
	2345
	2346	tauTES=0
	2347	do l=1,nlayermx
	2348	Qabsice = min(
	2349	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
	2350	& )
	2351	opTES(1,l)= 0.75 * Qabsice *
	2352	& zq(1,l,igcm_h2o_ice) *
[883]	2353	& (zplev(1,l) - zplev(1,l+1)) / g
[722]	2354	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
	2355	tauTES=tauTES+ opTES(1,l)
	2356	enddo
	2357	CALL WRITEDIAGFI(ngridmx,'tauTESap',
[358]	2358	& 'tau abs 825 cm-1',
[520]	2359	& '',0,tauTES)
[722]	2360	else
	2361
	2362	CALL WRITEDIAGFI(ngridmx,'tauTES',
[520]	2363	& 'tau abs 825 cm-1',
	2364	& '',0,taucloudtes)
[722]	2365	endif
[520]	2366
	2367	mtot = 0
	2368	icetot = 0
	2369	h2otot = qsurf(1,igcm_h2o_ice)
[722]	2370
[358]	2371	do l=1,nlayer
[520]	2372	mtot = mtot + zq(1,l,igcm_h2o_vap)
[883]	2373	& * (zplev(1,l) - zplev(1,l+1)) / g
[520]	2374	icetot = icetot + zq(1,l,igcm_h2o_ice)
[883]	2375	& * (zplev(1,l) - zplev(1,l+1)) / g
[722]	2376	end do
	2377	h2otot = h2otot+mtot+icetot
	2378
	2379	CALL WRITEDIAGFI(ngridmx,'h2otot',
	2380	& 'h2otot',
	2381	& 'kg/m2',0,h2otot)
	2382	CALL WRITEDIAGFI(ngridmx,'mtot',
	2383	& 'mtot',
	2384	& 'kg/m2',0,mtot)
	2385	CALL WRITEDIAGFI(ngridmx,'icetot',
	2386	& 'icetot',
	2387	& 'kg/m2',0,icetot)
	2388
	2389	if (scavenging) then
	2390
	2391	rave = 0
	2392	do l=1,nlayer
[669]	2393	cccc Column integrated effective ice radius
	2394	cccc is weighted by total ice surface area (BETTER)
[722]	2395	rave = rave + tauscaling(1) *
	2396	& zq(1,l,igcm_ccn_number) *
[883]	2397	& (zplev(1,l) - zplev(1,l+1)) / g *
[722]	2398	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
	2399	enddo
	2400	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
	2401
	2402	Nccntot= 0
[883]	2403	call watersat(ngridmx*nlayermx,zt,zplay,zqsat)
[520]	2404	do l=1,nlayermx
[669]	2405	Nccntot = Nccntot +
[520]	2406	& zq(1,l,igcm_ccn_number)*tauscaling(1)
[883]	2407	& *(zplev(1,l) - zplev(1,l+1)) / g
[520]	2408	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
	2409	satu(1,l) = (max(satu(1,l),float(1))-1)
	2410	! & * zq(1,l,igcm_h2o_vap) *
[883]	2411	! & (zplev(1,l) - zplev(1,l+1)) / g
[520]	2412	enddo
[722]	2413	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
	2414	& satu)
[669]	2415	CALL WRITEDIAGFI(ngridmx,'Nccntot',
	2416	& 'Nccntot',
	2417	& 'nbr/m2',0,Nccntot)
[411]	2418
[769]	2419	call WRITEDIAGFI(ngridmx,'zdqsed_dustq'
	2420	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
	2421	call WRITEDIAGFI(ngridmx,'zdqsed_dustN'
	2422	&,'sedimentation N','Nbr.m-2.s-1',1,
	2423	& zdqsed(1,:,igcm_dust_number))
[722]	2424
	2425	else ! of if (scavenging)
	2426
	2427	cccc Column integrated effective ice radius
	2428	cccc is weighted by total ice mass (LESS GOOD)
	2429	rave = 0
	2430	do l=1,nlayer
	2431	rave = rave + zq(1,l,igcm_h2o_ice)
[883]	2432	& * (zplev(1,l) - zplev(1,l+1)) / g
[722]	2433	& * rice(1,l) * (1.+nuice_ref)
	2434	enddo
	2435	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
	2436	endif ! of if (scavenging)
	2437
	2438	CALL WRITEDIAGFI(ngridmx,'reffice',
	2439	& 'reffice',
	2440	& 'm',0,rave)
	2441
	2442	do iq=1,nq
	2443	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_s',
	2444	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
	2445	end do
[358]	2446
[633]	2447	call WRITEDIAGFI(ngridmx,'zdqcloud_ice','cloud ice',
	2448	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
	2449	call WRITEDIAGFI(ngridmx,'zdqcloud_vap','cloud vap',
	2450	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
	2451	call WRITEDIAGFI(ngridmx,'zdqcloud','cloud ice',
	2452	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
	2453	& +zdqcloud(1,:,igcm_h2o_vap))
[411]	2454
[520]	2455	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
[411]	2456	& rice)
[520]	2457	ENDIF ! of IF (water)
[411]	2458
[358]	2459	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2460	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2461
	2462
[883]	2463	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
[226]	2464
	2465	do l=2,nlayer-1
	2466	tmean=zt(1,l)
	2467	if(zt(1,l).ne.zt(1,l-1))
	2468	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	2469	zlocal(l)= zlocal(l-1)
[883]	2470	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
[226]	2471	enddo
	2472	zlocal(nlayer)= zlocal(nlayer-1)-
[883]	2473	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
[226]	2474	& rnew(1,nlayer)*tmean/g
	2475
[42]	2476	END IF ! if(ngrid.ne.1)
	2477
	2478	icount=icount+1
	2479	RETURN
	2480	END

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

Download in other formats:

Original Format