Context Navigation

physiq.F @ 1242

Last change on this file since 1242 was 234, checked in by aslmd, 13 years ago

LMDZ.MARS
+ MESOSCALE

An important change to merge physiq.F and inifis.F for GCM and mesoscale
This is mostly transparent to GCM users and developers (use of MESOSCALE precompiler flags)
Makes it easy (and mostly automatic!) for changes in GCM physics to be impacted in mesoscale physics
A few minor changes have followed in the GCM (slope scheme, one-point diagnostic).

Compilation + run is OK on both sides (GCM and mesoscale).
On the mesoscale side, call_meso_physiq?.inc and call_meso_inifis?.inc have been changed accordingly.

Here is the excerpt from README file:

19/07/2011 == AS

Finished converging meso_physiq.F and meso_inifis.F towards physiq.F and inifis.F --> see previous point 15/07/2011 --> meso_ routines no longer exist (everything is in meso_inc and transparent to GCM users) --> GCM routines include mesoscale parts within MESOSCALE precompiler commands --> MESOSCALE parts are as hidden as possible not to mess up with GCM users/developers
Cleaned inelegant or useless #ifdef [or] #ifndef MESOSCALE in physiq and inifis so that a minimum amount of such precompiler commands is now reached [mainly related to I/O]
Added the SF08 slope insolation model in the general physics parameterizations. Added a callslope keyword in inifis.F and callkeys.h --> This keyword is False by default / True if you use -DMESOSCALE
Removed the obsolete call to Viking Lander 1 diagnostic Replaced it with a diagnostic for opacity at the domain center [valid for GCM and mesoscale]

File size: 61.6 KB

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

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

Download in other formats:

Original Format