Context Navigation

physiq.F @ 73

Last change on this file since 73 was 38, checked in by emillour, 15 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 52.6 KB

Rev	Line
[38]	1	SUBROUTINE physiq(ngrid,nlayer,nq,
	2	$ firstcall,lastcall,
	3	$ pday,ptime,ptimestep,
	4	$ pplev,pplay,pphi,
	5	$ pu,pv,pt,pq,
	6	$ pw,
	7	$ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
	8
	9
	10	IMPLICIT NONE
	11	c=======================================================================
	12	c
	13	c subject:
	14	c --------
	15	c
	16	c Organisation of the physical parametrisations of the LMD
	17	c martian atmospheric general circulation model.
	18	c
	19	c The GCM can be run without or with tracer transport
	20	c depending on the value of Logical "tracer" in file "callphys.def"
	21	c Tracers may be water vapor, ice OR chemical species OR dust particles
	22	c
	23	c SEE comments in initracer.F about numbering of tracer species...
	24	c
	25	c It includes:
	26	c
	27	c 1. Initialization:
	28	c 1.1 First call initializations
	29	c 1.2 Initialization for every call to physiq
	30	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	31	c 2. Compute radiative transfer tendencies
	32	c (longwave and shortwave) for CO2 and aerosols.
	33	c 3. Gravity wave and subgrid scale topography drag :
	34	c 4. Vertical diffusion (turbulent mixing):
	35	c 5. Convective adjustment
	36	c 6. Condensation and sublimation of carbon dioxide.
	37	c 7. TRACERS :
	38	c 7a. water and water ice
	39	c 7b. call for photochemistry when tracers are chemical species
	40	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
	41	c 7d. updates (CO2 pressure variations, surface budget)
	42	c 8. Contribution to tendencies due to thermosphere
	43	c 9. Surface and sub-surface temperature calculations
	44	c 10. Write outputs :
	45	c - "startfi", "histfi" (if it's time)
	46	c - Saving statistics (if "callstats = .true.")
	47	c - Dumping eof (if "calleofdump = .true.")
	48	c - Output any needed variables in "diagfi"
	49	c 11. Diagnostic: mass conservation of tracers
	50	c
	51	c author:
	52	c -------
	53	c Frederic Hourdin 15/10/93
	54	c Francois Forget 1994
	55	c Christophe Hourdin 02/1997
	56	c Subroutine completly rewritten by F.Forget (01/2000)
	57	c Introduction of the photochemical module: S. Lebonnois (11/2002)
	58	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
	59	c Water ice clouds: Franck Montmessin (update 06/2003)
	60	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	61	c Nb: See callradite.F for more information.
	62	c
	63	c arguments:
	64	c ----------
	65	c
	66	c input:
	67	c ------
	68	c ecri period (in dynamical timestep) to write output
	69	c ngrid Size of the horizontal grid.
	70	c All internal loops are performed on that grid.
	71	c nlayer Number of vertical layers.
	72	c nq Number of advected fields
	73	c firstcall True at the first call
	74	c lastcall True at the last call
	75	c pday Number of days counted from the North. Spring
	76	c equinoxe.
	77	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	78	c ptimestep timestep (s)
	79	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	80	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	81	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	82	c pu(ngrid,nlayer) u component of the wind (ms-1)
	83	c pv(ngrid,nlayer) v component of the wind (ms-1)
	84	c pt(ngrid,nlayer) Temperature (K)
	85	c pq(ngrid,nlayer,nq) Advected fields
	86	c pudyn(ngrid,nlayer) \
	87	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	88	c ptdyn(ngrid,nlayer) / corresponding variables
	89	c pqdyn(ngrid,nlayer,nq) /
	90	c pw(ngrid,?) vertical velocity
	91	c
	92	c output:
	93	c -------
	94	c
	95	c pdu(ngrid,nlayermx) \
	96	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
	97	c pdt(ngrid,nlayermx) / variables due to physical processes.
	98	c pdq(ngrid,nlayermx,nqmx) /
	99	c pdpsrf(ngrid) /
	100	c tracerdyn call tracer in dynamical part of GCM ?
	101
	102	c
	103	c=======================================================================
	104	c
	105	c 0. Declarations :
	106	c ------------------
	107
	108	#include "dimensions.h"
	109	#include "dimphys.h"
	110	#include "comgeomfi.h"
	111	#include "surfdat.h"
	112	#include "comsoil.h"
	113	#include "comdiurn.h"
	114	#include "callkeys.h"
	115	#include "comcstfi.h"
	116	#include "planete.h"
	117	#include "comsaison.h"
	118	#include "control.h"
	119	#include "dimradmars.h"
	120	#include "comg1d.h"
	121	#include "tracer.h"
	122	#include "nlteparams.h"
	123
	124	#include "chimiedata.h"
	125	#include "watercap.h"
	126	#include "param.h"
	127	#include "param_v3.h"
	128	#include "conc.h"
	129
	130	#include "netcdf.inc"
	131
	132
	133
	134	c Arguments :
	135	c -----------
	136
	137	c inputs:
	138	c -------
	139	INTEGER ngrid,nlayer,nq
	140	REAL ptimestep
	141	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	142	REAL pphi(ngridmx,nlayer)
	143	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	144	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	145	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
	146	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	147	LOGICAL firstcall,lastcall
	148
	149	REAL pday
	150	REAL ptime
	151	logical tracerdyn
	152
	153	c outputs:
	154	c --------
	155	c physical tendencies
	156	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	157	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	158	REAL pdpsrf(ngridmx) ! surface pressure tendency
	159
	160
	161	c Local saved variables:
	162	c ----------------------
	163	c aerosol (dust or ice) extinction optical depth at reference wavelength
	164	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
	165	c aerosol optical properties :
	166	REAL aerosol(ngridmx,nlayermx,naerkind)
	167
	168	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	169	INTEGER icount ! counter of calls to physiq during the run.
	170	REAL tsurf(ngridmx) ! Surface temperature (K)
	171	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	172	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
	173	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
	174	REAL emis(ngridmx) ! Thermal IR surface emissivity
	175	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	176	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	177	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	178	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	179	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
	180	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	181	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	182	INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic)
	183
	184	c Variables used by the water ice microphysical scheme:
	185	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
	186	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
	187	! of the size distribution
	188	c Albedo of deposited surface ice
	189	REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
	190
	191	SAVE day_ini, icount
	192	SAVE aerosol, tsurf,tsoil
	193	SAVE co2ice,albedo,emis, q2
	194	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
	195	SAVE ig_vl1
	196
	197	REAL stephan
	198	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	199	SAVE stephan
	200
	201	c Local variables :
	202	c -----------------
	203
	204	REAL CBRT
	205	EXTERNAL CBRT
	206
	207	CHARACTER*80 fichier
	208	INTEGER l,ig,ierr,igout,iq,i, tapphys
	209
	210	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
	211	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
	212	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
	213	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
	214	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
	215	! (used if active=F)
	216	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	217	REAL zls ! solar longitude (rad)
	218	REAL zday ! date (time since Ls=0, in martian days)
	219	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	220	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
	221	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
	222
	223	c Tendancies due to various processes:
	224	REAL dqsurf(ngridmx,nqmx)
	225	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	226	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
	227	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
	228	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
	229	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
	230	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
	231	REAL zdtsurf(ngridmx) ! (K/s)
	232	REAL zdtcloud(ngridmx,nlayermx)
	233	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	234	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	235	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	236	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	237	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	238	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	239	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	240	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	241
	242	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
	243	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
	244	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
	245	REAL zdqadj(ngridmx,nlayermx,nqmx)
	246	REAL zdqc(ngridmx,nlayermx,nqmx)
	247	REAL zdqcloud(ngridmx,nlayermx,nqmx)
	248	REAL zdqscloud(ngridmx,nqmx)
	249	REAL zdqchim(ngridmx,nlayermx,nqmx)
	250	REAL zdqschim(ngridmx,nqmx)
	251
	252	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	253	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	254	REAL zdumolvis(ngridmx,nlayermx)
	255	REAL zdvmolvis(ngridmx,nlayermx)
	256	real zdqmoldiff(ngridmx,nlayermx,nqmx)
	257
	258	c Local variable for local intermediate calcul:
	259	REAL zflubid(ngridmx)
	260	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	261	REAL zdum1(ngridmx,nlayermx)
	262	REAL zdum2(ngridmx,nlayermx)
	263	REAL ztim1,ztim2,ztim3, z1,z2
	264	REAL ztime_fin
	265	REAL zdh(ngridmx,nlayermx)
	266	INTEGER length
	267	PARAMETER (length=100)
	268
	269	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	270	c -----------------------------------------------------------------------------
	271	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	272	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	273	REAL zq(ngridmx,nlayermx,nqmx)
	274	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	275	character*2 str2
	276	character*5 str5
	277	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
	278	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
	279	! (particules kg-1)
	280	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
	281	real qtot1,qtot2 ! total aerosol mass
	282	integer igmin, lmin
	283	logical tdiag
	284
	285	real co2col(ngridmx) ! CO2 column
	286	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
	287	REAL zstress(ngrid), cd
	288	real hco2(nqmx),tmean, zlocal(nlayermx)
	289	real rho(ngridmx,nlayermx) ! density
	290	real vmr(ngridmx,nlayermx) ! volume mixing ratio
	291	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
	292	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
	293	REAL rave(ngridmx) ! Mean water ice effective radius (m)
	294	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
	295	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
	296	REAL Qabsice ! Water ice absorption coefficient
	297
	298
	299	REAL time_phys
	300
	301	c=======================================================================
	302
	303	c 1. Initialisation:
	304	c -----------------
	305
	306	c 1.1 Initialisation only at first call
	307	c ---------------------------------------
	308	IF (firstcall) THEN
	309
	310	c variables set to 0
	311	c ~~~~~~~~~~~~~~~~~~
	312	call zerophys(ngridnlayernaerkind,aerosol)
	313	call zerophys(ngrid*nlayer,dtrad)
	314	call zerophys(ngrid,fluxrad)
	315
	316	c read startfi
	317	c ~~~~~~~~~~~~
	318
	319	! Read netcdf initial physical parameters.
	320	CALL phyetat0 ("startfi.nc",0,0,
	321	& nsoilmx,nq,
	322	& day_ini,time_phys,
	323	& tsurf,tsoil,emis,q2,qsurf,co2ice)
	324
	325	if (pday.ne.day_ini) then
	326	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	327	& "physics and dynamics"
	328	write(,) "dynamics day: ",pday
	329	write(,) "physics day: ",day_ini
	330	stop
	331	endif
	332
	333	write (,) 'In physiq day_ini =', day_ini
	334
	335	c Initialize albedo and orbital calculation
	336	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	337	CALL surfini(ngrid,co2ice,qsurf,albedo)
	338	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	339
	340	c initialize soil
	341	c ~~~~~~~~~~~~~~~
	342	IF (callsoil) THEN
	343	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	344	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	345	ELSE
	346	PRINT*,
	347	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	348	DO ig=1,ngrid
	349	capcal(ig)=1.e5
	350	fluxgrd(ig)=0.
	351	ENDDO
	352	ENDIF
	353	icount=1
	354
	355
	356	c initialize tracers
	357	c ~~~~~~~~~~~~~~~~~~
	358	tracerdyn=tracer
	359	IF (tracer) THEN
	360	CALL initracer(qsurf,co2ice)
	361	ENDIF ! end tracer
	362
	363	c Determining gridpoint near Viking Lander 1 (used for diagnostic only)
	364	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	365
	366	if(ngrid.ne.1) then
	367	latvl1= 22.27
	368	lonvl1= -47.94
	369	ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim +
	370	& int(1.5+(lonvl1+180)*iim/360.)
	371	write(,) 'Viking Lander 1 GCM point: lat,lon',
	372	& lati(ig_vl1)180/pi, long(ig_vl1)180/pi
	373	end if
	374
	375	c Initialize thermospheric parameters
	376	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	377
	378	if (callthermos) call param_read
	379
	380	c Initialize R and Cp as constant
	381
	382	if (.not.callthermos .and. .not.photochem) then
	383	do l=1,nlayermx
	384	do ig=1,ngridmx
	385	rnew(ig,l)=r
	386	cpnew(ig,l)=cpp
	387	mmean(ig,l)=mugaz
	388	enddo
	389	enddo
	390	endif
	391
	392	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
	393	write(,)"physiq: water_param Surface ice alb:",alb_surfice
	394	ENDIF
	395
	396	ENDIF ! (end of "if firstcall")
	397
	398	c ---------------------------------------------------
	399	c 1.2 Initializations done at every physical timestep:
	400	c ---------------------------------------------------
	401	c
	402	IF (ngrid.NE.ngridmx) THEN
	403	PRINT*,'STOP in PHYSIQ'
	404	PRINT*,'Probleme de dimensions :'
	405	PRINT*,'ngrid = ',ngrid
	406	PRINT*,'ngridmx = ',ngridmx
	407	STOP
	408	ENDIF
	409
	410	c Initialize various variables
	411	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	412	call zerophys(ngrid*nlayer, pdv)
	413	call zerophys(ngrid*nlayer, pdu)
	414	call zerophys(ngrid*nlayer, pdt)
	415	call zerophys(ngridnlayernq, pdq)
	416	call zerophys(ngrid, pdpsrf)
	417	call zerophys(ngrid, zflubid)
	418	call zerophys(ngrid, zdtsurf)
	419	call zerophys(ngrid*nq, dqsurf)
	420	igout=ngrid/2+1
	421
	422
	423	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	424
	425	c Compute Solar Longitude (Ls) :
	426	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	427	if (season) then
	428	call solarlong(zday,zls)
	429	else
	430	call solarlong(float(day_ini),zls)
	431	end if
	432
	433	c Compute geopotential at interlayers
	434	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	435	c ponderation des altitudes au niveau des couches en dp/p
	436
	437	DO l=1,nlayer
	438	DO ig=1,ngrid
	439	zzlay(ig,l)=pphi(ig,l)/g
	440	ENDDO
	441	ENDDO
	442	DO ig=1,ngrid
	443	zzlev(ig,1)=0.
	444	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	445	ENDDO
	446	DO l=2,nlayer
	447	DO ig=1,ngrid
	448	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	449	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	450	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	451	ENDDO
	452	ENDDO
	453
	454
	455	! Potential temperature calculation not the same in physiq and dynamic
	456
	457	c Compute potential temperature
	458	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	459	DO l=1,nlayer
	460	DO ig=1,ngrid
	461	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	462	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	463	ENDDO
	464	ENDDO
	465
	466	c-----------------------------------------------------------------------
	467	c 1.2.5 Compute mean mass, cp, and R
	468	c --------------------------------
	469
	470	if(photochem.or.callthermos) then
	471	call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
	472	endif
	473	c-----------------------------------------------------------------------
	474	c 2. Compute radiative tendencies :
	475	c------------------------------------
	476
	477
	478	IF (callrad) THEN
	479	IF( MOD(icount-1,iradia).EQ.0) THEN
	480
	481	c Local Solar zenith angle
	482	c ~~~~~~~~~~~~~~~~~~~~~~~~
	483	CALL orbite(zls,dist_sol,declin)
	484
	485	IF(diurnal) THEN
	486	ztim1=SIN(declin)
	487	ztim2=COS(declin)COS(2.pi*(zday-.5))
	488	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	489
	490	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	491	s ztim1,ztim2,ztim3, mu0,fract)
	492
	493	ELSE
	494	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
	495	ENDIF
	496
	497	c NLTE cooling from CO2 emission
	498	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	499
	500	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
	501
	502	c Find number of layers for LTE radiation calculations
	503	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
	504	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
	505
	506	c Note: Dustopacity.F has been transferred to callradite.F
	507
	508	c Call main radiative transfer scheme
	509	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	510	c Transfer through CO2 (except NIR CO2 absorption)
	511	c and aerosols (dust and water ice)
	512
	513	c Radiative transfer
	514	c ------------------
	515	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
	516	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
	517	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
	518	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
	519
	520	c CO2 near infrared absorption
	521	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	522	call zerophys(ngrid*nlayer,zdtnirco2)
	523	if (callnirco2) then
	524	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
	525	. mu0,fract,declin, zdtnirco2)
	526	endif
	527
	528	c Radiative flux from the sky absorbed by the surface (W.m-2)
	529	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	530	DO ig=1,ngrid
	531	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	532	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
	533	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
	534	ENDDO
	535
	536
	537	c Net atmospheric radiative heating rate (K.s-1)
	538	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	539	IF(callnlte) THEN
	540	CALL blendrad(ngrid, nlayer, pplay,
	541	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
	542	ELSE
	543	DO l=1,nlayer
	544	DO ig=1,ngrid
	545	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	546	& +zdtnirco2(ig,l)
	547	ENDDO
	548	ENDDO
	549	ENDIF
	550
	551
	552
	553	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	554
	555	c Transformation of the radiative tendencies:
	556	c -------------------------------------------
	557
	558	c Net radiative surface flux (W.m-2)
	559	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	560	c
	561	DO ig=1,ngrid
	562	zplanck(ig)=tsurf(ig)*tsurf(ig)
	563	zplanck(ig)=emis(ig)*
	564	$ stephanzplanck(ig)zplanck(ig)
	565	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
	566	ENDDO
	567
	568
	569	DO l=1,nlayer
	570	DO ig=1,ngrid
	571	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	572	ENDDO
	573	ENDDO
	574
	575	ENDIF ! of IF (callrad)
	576
	577	c-----------------------------------------------------------------------
	578	c 3. Gravity wave and subgrid scale topography drag :
	579	c -------------------------------------------------
	580
	581
	582	IF(calllott)THEN
	583
	584	CALL calldrag_noro(ngrid,nlayer,ptimestep,
	585	& pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
	586
	587	DO l=1,nlayer
	588	DO ig=1,ngrid
	589	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
	590	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
	591	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
	592	ENDDO
	593	ENDDO
	594	ENDIF
	595
	596	c-----------------------------------------------------------------------
	597	c 4. Vertical diffusion (turbulent mixing):
	598	c -----------------------------------------
	599	c
	600	IF (calldifv) THEN
	601
	602
	603	DO ig=1,ngrid
	604	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	605	ENDDO
	606
	607	CALL zerophys(ngrid*nlayer,zdum1)
	608	CALL zerophys(ngrid*nlayer,zdum2)
	609	do l=1,nlayer
	610	do ig=1,ngrid
	611	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	612	enddo
	613	enddo
	614
	615	c Calling vdif (Martian version WITH CO2 condensation)
	616	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
	617	$ ptimestep,capcal,lwrite,
	618	$ pplay,pplev,zzlay,zzlev,z0,
	619	$ pu,pv,zh,pq,tsurf,emis,qsurf,
	620	$ zdum1,zdum2,zdh,pdq,zflubid,
	621	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
	622	& zdqdif,zdqsdif)
	623
	624	DO l=1,nlayer
	625	DO ig=1,ngrid
	626	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	627	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	628	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	629
	630	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	631
	632	ENDDO
	633	ENDDO
	634
	635	DO ig=1,ngrid
	636	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	637	ENDDO
	638
	639	if (tracer) then
	640	DO iq=1, nq
	641	DO l=1,nlayer
	642	DO ig=1,ngrid
	643	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	644	ENDDO
	645	ENDDO
	646	ENDDO
	647	DO iq=1, nq
	648	DO ig=1,ngrid
	649	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	650	ENDDO
	651	ENDDO
	652	end if ! of if (tracer)
	653
	654	ELSE
	655	DO ig=1,ngrid
	656	zdtsurf(ig)=zdtsurf(ig)+
	657	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	658	ENDDO
	659	ENDIF ! of IF (calldifv)
	660
	661
	662	c-----------------------------------------------------------------------
	663	c 5. Dry convective adjustment:
	664	c -----------------------------
	665
	666	IF(calladj) THEN
	667
	668	DO l=1,nlayer
	669	DO ig=1,ngrid
	670	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	671	ENDDO
	672	ENDDO
	673	CALL zerophys(ngrid*nlayer,zduadj)
	674	CALL zerophys(ngrid*nlayer,zdvadj)
	675	CALL zerophys(ngrid*nlayer,zdhadj)
	676
	677	CALL convadj(ngrid,nlayer,nq,ptimestep,
	678	$ pplay,pplev,zpopsk,
	679	$ pu,pv,zh,pq,
	680	$ pdu,pdv,zdh,pdq,
	681	$ zduadj,zdvadj,zdhadj,
	682	$ zdqadj)
	683
	684	DO l=1,nlayer
	685	DO ig=1,ngrid
	686	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	687	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	688	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	689
	690	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	691	ENDDO
	692	ENDDO
	693
	694	if(tracer) then
	695	DO iq=1, nq
	696	DO l=1,nlayer
	697	DO ig=1,ngrid
	698	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	699	ENDDO
	700	ENDDO
	701	ENDDO
	702	end if
	703	ENDIF ! of IF(calladj)
	704
	705	c-----------------------------------------------------------------------
	706	c 6. Carbon dioxide condensation-sublimation:
	707	c -------------------------------------------
	708
	709	IF (callcond) THEN
	710	CALL newcondens(ngrid,nlayer,nq,ptimestep,
	711	$ capcal,pplay,pplev,tsurf,pt,
	712	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	713	$ co2ice,albedo,emis,
	714	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
	715	$ fluxsurf_sw,zls)
	716
	717	DO l=1,nlayer
	718	DO ig=1,ngrid
	719	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	720	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	721	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	722	ENDDO
	723	ENDDO
	724	DO ig=1,ngrid
	725	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	726	ENDDO
	727
	728	IF (tracer) THEN
	729	DO iq=1, nq
	730	DO l=1,nlayer
	731	DO ig=1,ngrid
	732	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	733	ENDDO
	734	ENDDO
	735	ENDDO
	736	ENDIF ! of IF (tracer)
	737
	738	ENDIF ! of IF (callcond)
	739
	740	c-----------------------------------------------------------------------
	741	c 7. Specific parameterizations for tracers
	742	c: -----------------------------------------
	743
	744	if (tracer) then
	745
	746	c 7a. Water and ice
	747	c ---------------
	748
	749	c ---------------------------------------
	750	c Water ice condensation in the atmosphere
	751	c ----------------------------------------
	752	IF (water) THEN
	753
	754	call watercloud(ngrid,nlayer,ptimestep,
	755	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
	756	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
	757	& nq,naerkind,tau,
	758	& ccn,rdust,rice,nuice)
	759	if (activice) then
	760	c Temperature variation due to latent heat release
	761	DO l=1,nlayer
	762	DO ig=1,ngrid
	763	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
	764	ENDDO
	765	ENDDO
	766	endif
	767
	768	! increment water vapour and ice atmospheric tracers tendencies
	769	IF (water) THEN
	770	DO l=1,nlayer
	771	DO ig=1,ngrid
	772	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	773	& zdqcloud(ig,l,igcm_h2o_vap)
	774	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	775	& zdqcloud(ig,l,igcm_h2o_ice)
	776	ENDDO
	777	ENDDO
	778	ENDIF ! of IF (water) THEN
	779	! Increment water ice surface tracer tendency
	780	DO ig=1,ngrid
	781	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
	782	& zdqscloud(ig,igcm_h2o_ice)
	783	ENDDO
	784
	785	END IF ! of IF (water)
	786
	787	c 7b. Chemical species
	788	c ------------------
	789
	790	c --------------
	791	c photochemistry :
	792	c --------------
	793	IF (photochem .or. thermochem) then
	794	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
	795	& zzlay,zday,pq,pdq,rice,
	796	& zdqchim,zdqschim,zdqcloud,zdqscloud)
	797	!NB: Photochemistry includes condensation of H2O2
	798
	799	! increment values of tracers:
	800	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	801	! tracers is zero anyways
	802	DO l=1,nlayer
	803	DO ig=1,ngrid
	804	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	805	ENDDO
	806	ENDDO
	807	ENDDO ! of DO iq=1,nq
	808	! add condensation tendency for H2O2
	809	if (igcm_h2o2.ne.0) then
	810	DO l=1,nlayer
	811	DO ig=1,ngrid
	812	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	813	& +zdqcloud(ig,l,igcm_h2o2)
	814	ENDDO
	815	ENDDO
	816	endif
	817
	818	! increment surface values of tracers:
	819	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	820	! tracers is zero anyways
	821	DO ig=1,ngrid
	822	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	823	ENDDO
	824	ENDDO ! of DO iq=1,nq
	825	! add condensation tendency for H2O2
	826	if (igcm_h2o2.ne.0) then
	827	DO ig=1,ngrid
	828	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	829	& +zdqscloud(ig,igcm_h2o2)
	830	ENDDO
	831	endif
	832
	833	END IF ! of IF (photochem.or.thermochem)
	834
	835	c 7c. Aerosol particles
	836	c -------------------
	837
	838	c ----------
	839	c Dust devil :
	840	c ----------
	841	IF(callddevil) then
	842	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
	843	& zdqdev,zdqsdev)
	844
	845	if (dustbin.ge.1) then
	846	do iq=1,nq
	847	DO l=1,nlayer
	848	DO ig=1,ngrid
	849	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(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) + zdqsdev(ig,iq)
	856	ENDDO
	857	enddo
	858	endif ! of if (dustbin.ge.1)
	859
	860	END IF ! of IF (callddevil)
	861
	862	c -------------
	863	c Sedimentation : acts also on water ice
	864	c -------------
	865	IF (sedimentation) THEN
	866	!call zerophys(ngridnlayernq, zdqsed)
	867	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	868	!call zerophys(ngrid*nq, zdqssed)
	869	zdqssed(1:ngrid,1:nq)=0
	870
	871	call callsedim(ngrid,nlayer, ptimestep,
	872	& pplev,zzlev, pt, rdust, rice,
	873	& pq, pdq, zdqsed, zdqssed,nq)
	874
	875	DO iq=1, nq
	876	DO l=1,nlayer
	877	DO ig=1,ngrid
	878	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	879	ENDDO
	880	ENDDO
	881	ENDDO
	882	DO iq=1, nq
	883	DO ig=1,ngrid
	884	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	885	ENDDO
	886	ENDDO
	887	END IF ! of IF (sedimentation)
	888
	889	c 7d. Updates
	890	c ---------
	891
	892	DO iq=1, nq
	893	DO ig=1,ngrid
	894
	895	c ---------------------------------
	896	c Updating tracer budget on surface
	897	c ---------------------------------
	898	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	899
	900	ENDDO ! (ig)
	901	ENDDO ! (iq)
	902
	903	endif ! of if (tracer)
	904
	905
	906	c-----------------------------------------------------------------------
	907	c 8. THERMOSPHERE CALCULATION
	908	c-----------------------------------------------------------------------
	909
	910	if (callthermos) then
	911	call thermosphere(pplev,pplay,dist_sol,
	912	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	913	& pt,pq,pu,pv,pdt,pdq,
	914	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	915
	916	DO l=1,nlayer
	917	DO ig=1,ngrid
	918	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	919	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	920	& +zdteuv(ig,l)
	921	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	922	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	923	DO iq=1, nq
	924	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	925	ENDDO
	926	ENDDO
	927	ENDDO
	928
	929	endif ! of if (callthermos)
	930
	931	c-----------------------------------------------------------------------
	932	c 9. Surface and sub-surface soil temperature
	933	c-----------------------------------------------------------------------
	934	c
	935	c
	936	c 9.1 Increment Surface temperature:
	937	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	938
	939	DO ig=1,ngrid
	940	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	941	ENDDO
	942
	943	c Prescribe a cold trap at south pole (except at high obliquity !!)
	944	c Temperature at the surface is set there to be the temperature
	945	c corresponding to equilibrium temperature between phases of CO2
	946
	947	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
	948	if (caps.and.(obliquit.lt.27.)) then
	949	! NB: Updated surface pressure, at grid point 'ngrid', is
	950	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
	951	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	952	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	953	endif
	954	c -------------------------------------------------------------
	955	c Change of surface albedo (set to 0.4) in case of ground frost
	956	c everywhere except on the north permanent cap and in regions
	957	c covered by dry ice.
	958	c ALWAYS PLACE these lines after newcondens !!!
	959	c -------------------------------------------------------------
	960	do ig=1,ngrid
	961	if ((co2ice(ig).eq.0).and.
	962	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
	963	albedo(ig,1) = alb_surfice
	964	albedo(ig,2) = alb_surfice
	965	endif
	966	enddo ! of do ig=1,ngrid
	967	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
	968
	969	c
	970	c 9.2 Compute soil temperatures and subsurface heat flux:
	971	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	972	IF (callsoil) THEN
	973	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
	974	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	975	ENDIF
	976
	977	c-----------------------------------------------------------------------
	978	c 10. Write output files
	979	c ----------------------
	980
	981	c -------------------------------
	982	c Dynamical fields incrementation
	983	c -------------------------------
	984	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	985	! temperature, zonal and meridional wind
	986	DO l=1,nlayer
	987	DO ig=1,ngrid
	988	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	989	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	990	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	991	ENDDO
	992	ENDDO
	993
	994	! tracers
	995	DO iq=1, nq
	996	DO l=1,nlayer
	997	DO ig=1,ngrid
	998	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	999	ENDDO
	1000	ENDDO
	1001	ENDDO
	1002
	1003	! surface pressure
	1004	DO ig=1,ngrid
	1005	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1006	ENDDO
	1007
	1008	! pressure
	1009	DO l=1,nlayer
	1010	DO ig=1,ngrid
	1011	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1012	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1013	ENDDO
	1014	ENDDO
	1015
	1016	! Density
	1017	DO l=1,nlayer
	1018	DO ig=1,ngrid
	1019	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1020	ENDDO
	1021	ENDDO
	1022
	1023	c Compute surface stress : (NB: z0 is a common in planete.h)
	1024	c DO ig=1,ngrid
	1025	c cd = (0.4/log(zzlay(ig,1)/z0))**2
	1026	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1027	c ENDDO
	1028
	1029	c Sum of fluxes in solar spectral bands (for output only)
	1030	DO ig=1,ngrid
	1031	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1032	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
	1033	ENDDO
	1034	c ***** TEST ****************************************************
	1035	ztim1 = 999
	1036	DO l=1,nlayer
	1037	DO ig=1,ngrid
	1038	if (pt(ig,l).lt.ztim1) then
	1039	ztim1 = pt(ig,l)
	1040	igmin = ig
	1041	lmin = l
	1042	end if
	1043	ENDDO
	1044	ENDDO
	1045	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
	1046	write(,) 'PHYSIQ: stability WARNING :'
	1047	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1048	& 'ig l =', igmin, lmin
	1049	end if
	1050	c *******************************************************************
	1051
	1052	c ---------------------
	1053	c Outputs to the screen
	1054	c ---------------------
	1055
	1056	IF (lwrite) THEN
	1057	PRINT*,'Global diagnostics for the physics'
	1058	PRINT,'Variables and their increments x and dx/dt dt'
	1059	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1060	WRITE(*,'(2f10.5,2f15.5)')
	1061	s tsurf(igout),zdtsurf(igout)*ptimestep,
	1062	s pplev(igout,1),pdpsrf(igout)*ptimestep
	1063	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1064	WRITE(*,'(i4,6f10.5)') (l,
	1065	s pu(igout,l),pdu(igout,l)*ptimestep,
	1066	s pv(igout,l),pdv(igout,l)*ptimestep,
	1067	s pt(igout,l),pdt(igout,l)*ptimestep,
	1068	s l=1,nlayer)
	1069	ENDIF ! of IF (lwrite)
	1070
	1071	IF (ngrid.NE.1) THEN
	1072	print,'Ls =',zls180./pi,
	1073	& ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
	1074	& ' tau(Viking1) =',tau(ig_vl1,1)
	1075
	1076
	1077	c -------------------------------------------------------------------
	1078	c Writing NetCDF file "RESTARTFI" at the end of the run
	1079	c -------------------------------------------------------------------
	1080	c Note: 'restartfi' is stored just before dynamics are stored
	1081	c in 'restart'. Between now and the writting of 'restart',
	1082	c there will have been the itau=itau+1 instruction and
	1083	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1084	c thus we store for time=time+dtvr
	1085
	1086	IF(lastcall) THEN
	1087	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1088	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1089	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1090	. ptimestep,pday,
	1091	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
	1092	. area,albedodat,inertiedat,zmea,zstd,zsig,
	1093	. zgam,zthe)
	1094	ENDIF
	1095
	1096	c -------------------------------------------------------------------
	1097	c Calculation of diagnostic variables written in both stats and
	1098	c diagfi files
	1099	c -------------------------------------------------------------------
	1100
	1101	if (tracer) then
	1102	if (water) then
	1103
	1104	call zerophys(ngrid,mtot)
	1105	call zerophys(ngrid,icetot)
	1106	call zerophys(ngrid,rave)
	1107	call zerophys(ngrid,tauTES)
	1108	do ig=1,ngrid
	1109	do l=1,nlayermx
	1110	mtot(ig) = mtot(ig) +
	1111	& zq(ig,l,igcm_h2o_vap) *
	1112	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1113	icetot(ig) = icetot(ig) +
	1114	& zq(ig,l,igcm_h2o_ice) *
	1115	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1116	rave(ig) = rave(ig) +
	1117	& zq(ig,l,igcm_h2o_ice) *
	1118	& (pplev(ig,l) - pplev(ig,l+1)) / g *
	1119	& rice(ig,l) * (1.+nuice_ref)
	1120	c Computing abs optical depth at 825 cm-1 in each
	1121	c layer to simulate NEW TES retrieval
	1122	Qabsice = min(
	1123	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1124	& )
	1125	opTES(ig,l)= 0.75 * Qabsice *
	1126	& zq(ig,l,igcm_h2o_ice) *
	1127	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1128	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1129	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1130	enddo
	1131	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
	1132	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1133	enddo
	1134
	1135	endif ! of if (water)
	1136	endif ! of if (tracer)
	1137
	1138	c -----------------------------------------------------------------
	1139	c WSTATS: Saving statistics
	1140	c -----------------------------------------------------------------
	1141	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1142	c which can later be used to make the statistic files of the run:
	1143	c "stats") only possible in 3D runs !
	1144
	1145
	1146	IF (callstats) THEN
	1147
	1148	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1149	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1150	call wstats(ngrid,"co2ice","CO2 ice cover",
	1151	& "kg.m-2",2,co2ice)
	1152	call wstats(ngrid,"fluxsurf_lw",
	1153	& "Thermal IR radiative flux to surface","W.m-2",2,
	1154	& fluxsurf_lw)
	1155	call wstats(ngrid,"fluxsurf_sw",
	1156	& "Solar radiative flux to surface","W.m-2",2,
	1157	& fluxsurf_sw_tot)
	1158	call wstats(ngrid,"fluxtop_lw",
	1159	& "Thermal IR radiative flux to space","W.m-2",2,
	1160	& fluxtop_lw)
	1161	call wstats(ngrid,"fluxtop_sw",
	1162	& "Solar radiative flux to space","W.m-2",2,
	1163	& fluxtop_sw_tot)
	1164	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1165	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1166	call wstats(ngrid,"v","Meridional (North-South) wind",
	1167	& "m.s-1",3,zv)
	1168	call wstats(ngrid,"w","Vertical (down-up) wind",
	1169	& "m.s-1",3,pw)
	1170	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
	1171	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
	1172	c call wstats(ngrid,"q2",
	1173	c & "Boundary layer eddy kinetic energy",
	1174	c & "m2.s-2",3,q2)
	1175	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1176	c & emis)
	1177	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1178	c & 2,zstress)
	1179	c call wstats(ngrid,"sw_htrt","sw heat.rate",
	1180	c & "W.m-2",3,zdtsw)
	1181	c call wstats(ngrid,"lw_htrt","lw heat.rate",
	1182	c & "W.m-2",3,zdtlw)
	1183
	1184	if (tracer) then
	1185	if (water) then
	1186	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1187	& *mugaz/mmol(igcm_h2o_vap)
	1188	call wstats(ngrid,"vmr_h2ovapor",
	1189	& "H2O vapor volume mixing ratio","mol/mol",
	1190	& 3,vmr)
	1191	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1192	& *mugaz/mmol(igcm_h2o_ice)
	1193	call wstats(ngrid,"vmr_h2oice",
	1194	& "H2O ice volume mixing ratio","mol/mol",
	1195	& 3,vmr)
	1196	call wstats(ngrid,"h2o_ice_s",
	1197	& "surface h2o_ice","kg/m2",
	1198	& 2,qsurf(1,igcm_h2o_ice))
	1199
	1200	call wstats(ngrid,"mtot",
	1201	& "total mass of water vapor","kg/m2",
	1202	& 2,mtot)
	1203	call wstats(ngrid,"icetot",
	1204	& "total mass of water ice","kg/m2",
	1205	& 2,icetot)
	1206	call wstats(ngrid,"reffice",
	1207	& "Mean reff","m",
	1208	& 2,rave)
	1209	c call wstats(ngrid,"rice",
	1210	c & "Ice particle size","m",
	1211	c & 3,rice)
	1212	c If activice is true, tauTES is computed in aeropacity.F;
	1213	if (.not.activice) then
	1214	call wstats(ngrid,"tauTESap",
	1215	& "tau abs 825 cm-1","",
	1216	& 2,tauTES)
	1217	endif
	1218
	1219	endif ! of if (water)
	1220
	1221	if (thermochem.or.photochem) then
	1222	do iq=1,nq
	1223	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
	1224	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
	1225	. (noms(iq).eq."h2").or.
	1226	. (noms(iq).eq."o3")) then
	1227	do l=1,nlayer
	1228	do ig=1,ngrid
	1229	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
	1230	end do
	1231	end do
	1232	call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1233	. "Volume mixing ratio","mol/mol",3,vmr)
	1234	endif
	1235	enddo
	1236	endif ! of if (thermochem.or.photochem)
	1237
	1238	endif ! of if (tracer)
	1239
	1240	IF(lastcall) THEN
	1241	write (,) "Writing stats..."
	1242	call mkstats(ierr)
	1243	ENDIF
	1244
	1245	ENDIF !if callstats
	1246
	1247	c (Store EOF for Mars Climate database software)
	1248	IF (calleofdump) THEN
	1249	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1250	ENDIF
	1251
	1252	c ==========================================================
	1253	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1254	c any variable for diagnostic (output with period
	1255	c "ecritphy", set in "run.def")
	1256	c ==========================================================
	1257	c WRITEDIAGFI can ALSO be called from any other subroutines
	1258	c for any variables !!
	1259	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1260	c & emis)
	1261	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1262	& tsurf)
	1263	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1264	c call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1265	c & co2ice)
	1266	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1267	c & "K",2,zt(1,7))
	1268	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1269	& fluxsurf_lw)
	1270	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1271	& fluxsurf_sw_tot)
	1272	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1273	& fluxtop_lw)
	1274	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1275	& fluxtop_sw_tot)
	1276	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1277	c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1278	c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1279	c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1280	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
	1281	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
	1282	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
	1283	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
	1284	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1285	c & zstress)
	1286	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	1287	c & 'w.m-2',3,zdtsw)
	1288	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	1289	c & 'w.m-2',3,zdtlw)
	1290
	1291	c ----------------------------------------------------------
	1292	c Outputs of the CO2 cycle
	1293	c ----------------------------------------------------------
	1294
	1295	if (tracer.and.(igcm_co2.ne.0)) then
	1296	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	1297	! & "kg/kg",2,zq(1,1,igcm_co2))
	1298	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	1299	! & "kg/kg",3,zq(1,1,igcm_co2))
	1300
	1301	! Compute co2 column
	1302	call zerophys(ngrid,co2col)
	1303	do l=1,nlayermx
	1304	do ig=1,ngrid
	1305	co2col(ig)=co2col(ig)+
	1306	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
	1307	enddo
	1308	enddo
	1309	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	1310	& co2col)
	1311	endif ! of if (tracer.and.(igcm_co2.ne.0))
	1312
	1313	c ----------------------------------------------------------
	1314	c Outputs of the water cycle
	1315	c ----------------------------------------------------------
	1316	if (tracer) then
	1317	if (water) then
	1318
	1319	CALL WRITEDIAGFI(ngridmx,'mtot',
	1320	& 'total mass of water vapor',
	1321	& 'kg/m2',2,mtot)
	1322	CALL WRITEDIAGFI(ngridmx,'icetot',
	1323	& 'total mass of water ice',
	1324	& 'kg/m2',2,icetot)
	1325	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1326	c & *mugaz/mmol(igcm_h2o_ice)
	1327	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
	1328	c & 'mol/mol',3,vmr)
	1329	CALL WRITEDIAGFI(ngridmx,'reffice',
	1330	& 'Mean reff',
	1331	& 'm',2,rave)
	1332	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
	1333	c & 'm',3,rice)
	1334	c If activice is true, tauTES is computed in aeropacity.F;
	1335	if (.not.activice) then
	1336	CALL WRITEDIAGFI(ngridmx,'tauTESap',
	1337	& 'tau abs 825 cm-1',
	1338	& '',2,tauTES)
	1339	endif
	1340	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
	1341	& 'surface h2o_ice',
	1342	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
	1343	endif !(water)
	1344
	1345
	1346	if (water.and..not.photochem) then
	1347	iq=nq
	1348	c write(str2(1:2),'(i2.2)') iq
	1349	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
	1350	c & 'kg.m-2',2,zdqscloud(1,iq))
	1351	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
	1352	c & 'kg/kg',3,zdqchim(1,1,iq))
	1353	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
	1354	c & 'kg/kg',3,zdqdif(1,1,iq))
	1355	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
	1356	c & 'kg/kg',3,zdqadj(1,1,iq))
	1357	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
	1358	c & 'kg/kg',3,zdqc(1,1,iq))
	1359	endif !(water.and..not.photochem)
	1360	endif
	1361
	1362	c ----------------------------------------------------------
	1363	c Outputs of the dust cycle
	1364	c ----------------------------------------------------------
	1365
	1366	c call WRITEDIAGFI(ngridmx,'tauref',
	1367	c & 'Dust ref opt depth','NU',2,tauref)
	1368
	1369	if (tracer.and.(dustbin.ne.0)) then
	1370	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
	1371	if (doubleq) then
	1372	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
	1373	c & 'kg.m-2',2,qsurf(1,1))
	1374	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
	1375	c & 'N.m-2',2,qsurf(1,2))
	1376	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
	1377	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
	1378	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
	1379	c & 'kg.m-2.s-1',2,zdqssed(1,1))
	1380	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
	1381	& 'm',3,rdust*ref_r0)
	1382	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
	1383	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
	1384	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
	1385	c & 'part/kg',3,pq(1,1,igcm_dust_number))
	1386	else
	1387	do iq=1,dustbin
	1388	write(str2(1:2),'(i2.2)') iq
	1389	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
	1390	& 'kg/kg',3,zq(1,1,iq))
	1391	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
	1392	& 'kg.m-2',2,qsurf(1,iq))
	1393	end do
	1394	endif ! (doubleq)
	1395	c if (submicron) then
	1396	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
	1397	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
	1398	c endif ! (submicron)
	1399	end if ! (tracer.and.(dustbin.ne.0))
	1400
	1401	c ----------------------------------------------------------
	1402	c Output in netcdf file "diagsoil.nc" for subterranean
	1403	c variables (output every "ecritphy", as for writediagfi)
	1404	c ----------------------------------------------------------
	1405
	1406	! Write soil temperature
	1407	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	1408	! & 3,tsoil)
	1409	! Write surface temperature
	1410	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	1411	! & 2,tsurf)
	1412
	1413	c ==========================================================
	1414	c END OF WRITEDIAGFI
	1415	c ==========================================================
	1416
	1417	ELSE ! if(ngrid.eq.1)
	1418
	1419	print,'Ls =',zls180./pi,
	1420	& ' tauref(700 Pa) =',tauref
	1421	c ----------------------------------------------------------------------
	1422	c Output in grads file "g1d" (ONLY when using testphys1d)
	1423	c (output at every X physical timestep)
	1424	c ----------------------------------------------------------------------
	1425	c
	1426	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	1427	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	1428	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	1429
	1430	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	1431	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	1432	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	1433	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	1434
	1435	! or output in diagfi.nc (for testphys1d)
	1436	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
	1437	call WRITEDIAGFI(ngridmx,'temp','Temperature',
	1438	& 'K',1,zt)
	1439
	1440	if(tracer) then
	1441	c CALL writeg1d(ngrid,1,tau,'tau','SI')
	1442	do iq=1,nq
	1443	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	1444	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
	1445	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	1446	end do
	1447	end if
	1448
	1449	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
	1450
	1451	do l=2,nlayer-1
	1452	tmean=zt(1,l)
	1453	if(zt(1,l).ne.zt(1,l-1))
	1454	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	1455	zlocal(l)= zlocal(l-1)
	1456	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
	1457	enddo
	1458	zlocal(nlayer)= zlocal(nlayer-1)-
	1459	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
	1460	& rnew(1,nlayer)*tmean/g
	1461
	1462	END IF ! if(ngrid.ne.1)
	1463
	1464	icount=icount+1
	1465	RETURN
	1466	END

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

Download in other formats:

Original Format