Context Navigation

physiq.F @ 1268

Last change on this file since 1268 was 1268, checked in by aslmd, 11 years ago
LMDZ.MARS. related to r1266. forgot to remove a few now-obsolete dimensions.h includes in Mars physics.
File size: 97.6 KB

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

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

Download in other formats:

Original Format