Context Navigation

physiq.F @ 1467

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

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

Download in other formats:

Original Format