Context Navigation

physiq.F @ 1307

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

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

Download in other formats:

Original Format