Context Navigation

physiq.F @ 135

Last change on this file since 135 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 47.7 KB

Rev	Line
[135]	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	use radinc_h, only : naerkind
	10
	11	implicit none
	12
	13
	14	!==================================================================
	15	!
	16	! Purpose
	17	! -------
	18	! Central subroutine for all the physics parameterisations in the
	19	! universal model. Originally adapted from the Mars LMDZ model.
	20	!
	21	! The model can be run without or with tracer transport
	22	! depending on the value of "tracer" in file "callphys.def".
	23	!
	24	!
	25	! It includes:
	26	!
	27	! 1. Initialization:
	28	! 1.1 Firstcall initializations
	29	! 1.2 Initialization for every call to physiq
	30	! 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	31	! 2. Compute radiative transfer tendencies
	32	! (longwave and shortwave).
	33	! 4. Vertical diffusion (turbulent mixing):
	34	! 5. Convective adjustment
	35	! 6. Condensation and sublimation of gases (currently just CO2).
	36	! 7. TRACERS :
	37	! 7a. water and water ice
	38	! 7c. other schemes for tracer transport (lifting, sedimentation)
	39	! 7d. updates (pressure variations, surface budget)
	40	! 9. Surface and sub-surface temperature calculations
	41	! 10. Write outputs :
	42	! - "startfi", "histfi" if it's time
	43	! - Saving statistics if "callstats = .true."
	44	! - Output any needed variables in "diagfi"
	45	! 10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
	46	!
	47	! arguments
	48	! ---------
	49	!
	50	! input
	51	! -----
	52	! ecri period (in dynamical timestep) to write output
	53	! ngrid Size of the horizontal grid.
	54	! All internal loops are performed on that grid.
	55	! nlayer Number of vertical layers.
	56	! nq Number of advected fields
	57	! firstcall True at the first call
	58	! lastcall True at the last call
	59	! pday Number of days counted from the North. Spring
	60	! equinoxe.
	61	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	62	! ptimestep timestep (s)
	63	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	64	! pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	65	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	66	! pu(ngrid,nlayer) u component of the wind (ms-1)
	67	! pv(ngrid,nlayer) v component of the wind (ms-1)
	68	! pt(ngrid,nlayer) Temperature (K)
	69	! pq(ngrid,nlayer,nq) Advected fields
	70	! pudyn(ngrid,nlayer) \
	71	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	72	! ptdyn(ngrid,nlayer) / corresponding variables
	73	! pqdyn(ngrid,nlayer,nq) /
	74	! pw(ngrid,?) vertical velocity
	75	!
	76	! output
	77	! ------
	78	!
	79	! pdu(ngrid,nlayermx) \
	80	! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
	81	! pdt(ngrid,nlayermx) / variables due to physical processes.
	82	! pdq(ngrid,nlayermx) /
	83	! pdpsrf(ngrid) /
	84	! tracerdyn call tracer in dynamical part of GCM ?
	85	!
	86	!
	87	! Authors
	88	! -------
	89	! Frederic Hourdin 15/10/93
	90	! Francois Forget 1994
	91	! Christophe Hourdin 02/1997
	92	! Subroutine completly rewritten by F.Forget (01/2000)
	93	! Water ice clouds: Franck Montmessin (update 06/2003)
	94	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	95	! New correlated-k radiative scheme: R. Wordsworth (2009)
	96	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
	97	!
	98	!==================================================================
	99
	100
	101	c 0. Declarations :
	102	c ------------------
	103
	104	#include "dimensions.h"
	105	#include "dimphys.h"
	106	#include "comgeomfi.h"
	107	#include "surfdat.h"
	108	#include "comsoil.h"
	109	#include "comdiurn.h"
	110	#include "callkeys.h"
	111	#include "comcstfi.h"
	112	#include "planete.h"
	113	#include "comsaison.h"
	114	#include "control.h"
	115	#include "comg1d.h"
	116	#include "tracer.h"
	117
	118	#include "watercap.h"
	119
	120	#include "netcdf.inc"
	121
	122
	123
	124	c Arguments :
	125	c -----------
	126
	127	c inputs:
	128	c -------
	129	INTEGER ngrid,nlayer,nq
	130	REAL ptimestep
	131	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	132	REAL pphi(ngridmx,nlayer)
	133	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	134	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	135	REAL pw(ngridmx,nlayer) ! pvervel transmitted by dyn3d
	136	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	137	LOGICAL firstcall,lastcall
	138
	139	REAL pday
	140	REAL ptime
	141	logical tracerdyn
	142
	143	c outputs:
	144	c --------
	145	c physical tendencies
	146	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	147	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	148	REAL pdpsrf(ngridmx) ! surface pressure tendency
	149
	150
	151	c Local saved variables:
	152	c ----------------------
	153	c aerosol (dust or ice) extinction optical depth at reference wavelength
	154	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
	155	c aerosol optical properties:
	156	c REAL aerosol(ngridmx,nlayermx,naerkind)
	157	! this is now internal to callcorrk and hence no longer needed here
	158
	159	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	160	INTEGER icount ! counter of calls to physiq during the run.
	161	REAL tsurf(ngridmx) ! Surface temperature (K)
	162	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	163	REAL albedo(ngridmx) ! Surface albedo
	164
	165	real albedo0(ngridmx) ! Surface albedo
	166	save albedo0
	167
	168	REAL emis(ngridmx) ! Thermal IR surface emissivity
	169	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	170	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	171	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	172	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	173	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
	174	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	175	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	176
	177	SAVE day_ini, icount
	178	SAVE tsurf,tsoil
	179	SAVE albedo,emis, q2
	180	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
	181
	182	REAL stephan
	183	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	184	SAVE stephan
	185
	186	c Local variables :
	187	c -----------------
	188
	189	! aerosol (dust or ice) extinction optical depth at reference wavelength
	190	! for the "naerkind" optically active aerosols:
	191	REAL aerosol(ngridmx,nlayermx,naerkind)
	192
	193	CHARACTER*80 fichier
	194	INTEGER l,ig,ierr,igout,iq,i, tapphys
	195	! INTEGER iqmin ! Used if iceparty engaged
	196
	197	REAL fluxsurf_lw(ngridmx) ! incident LW (IR) surface flux (W.m-2)
	198	REAL fluxsurf_sw(ngridmx) ! incident SW (solar) surface flux (W.m-2)
	199	REAL fluxtop_lw(ngridmx) ! Outgoing LW (IR) flux to space (W.m-2)
	200	REAL fluxtop_sw(ngridmx) ! Outgoing SW (solar) flux to space (W.m-2)
	201
	202	! included by RW for equilibration test
	203	real fluxtop_dn(ngridmx)
	204	real fluxabs_sw(ngridmx) ! absorbed shortwave radiation
	205	real fluxdyn(ngridmx) ! horizontal heat transport by dynamics
	206
	207	REAL zls ! solar longitude (rad)
	208	REAL zday ! date (time since Ls=0, in martian days)
	209	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	210	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
	211	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
	212
	213	c Aerosol effective radius used for radiative transfer (units=meters)
	214	REAL reffrad(ngridmx,nlayermx,naerkind)
	215
	216	c Tendencies due to various processes:
	217	REAL dqsurf(ngridmx,nqmx)
	218	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	219	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
	220	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear areas
	221	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear areas
	222	REAL zdtsurf(ngridmx) ! (K/s)
	223	REAL zdtlscale(ngridmx,nlayermx)
	224	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	225	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	226	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	227	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	228	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	229	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	230	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	231	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	232
	233	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
	234	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
	235	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
	236	REAL zdqadj(ngridmx,nlayermx,nqmx)
	237	REAL zdqc(ngridmx,nlayermx,nqmx)
	238	REAL zdqlscale(ngridmx,nlayermx,nqmx)
	239	REAL zdqslscale(ngridmx,nqmx)
	240	REAL zdqchim(ngridmx,nlayermx,nqmx)
	241	REAL zdqschim(ngridmx,nqmx)
	242
	243	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	244	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	245	REAL zdumolvis(ngridmx,nlayermx)
	246	REAL zdvmolvis(ngridmx,nlayermx)
	247	real zdqmoldiff(ngridmx,nlayermx,nqmx)
	248
	249	c Local variable for local intermediate calcul:
	250	REAL zflubid(ngridmx)
	251	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	252	REAL zdum1(ngridmx,nlayermx)
	253	REAL zdum2(ngridmx,nlayermx)
	254	REAL ztim1,ztim2,ztim3, z1,z2
	255	REAL ztime_fin
	256	REAL zdh(ngridmx,nlayermx)
	257	INTEGER length
	258	PARAMETER (length=100)
	259
	260	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	261	c -----------------------------------------------------------------------------
	262	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	263	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	264	REAL zq(ngridmx,nlayermx,nqmx)
	265	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	266	character*2 str2
	267	character*5 str5
	268	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
	269	real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
	270	save ztprevious
	271	real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
	272	real qtot1,qtot2 ! total aerosol mass
	273	integer igmin, lmin
	274	logical tdiag
	275
	276	real co2col(ngridmx) ! CO2 column
	277	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
	278	REAL zstress(ngrid), cd
	279	real hco2(nqmx),tmean, zlocal(nlayermx)
	280	real rho(ngridmx,nlayermx) ! density
	281	real vmr(ngridmx,nlayermx) ! volume mixing ratio
	282
	283	REAL time_phys
	284
	285	! included by RW for kastprof scheme
	286	logical kastprof
	287
	288	! reinstated by RW for diagnostic
	289	REAL tau_col(ngrid)
	290
	291	! included by RW for tidally locked dynamics
	292	! logical tlocked
	293
	294	! included by RW to reduce insanity of code
	295	real flux1,flux2,flux3,ts1,ts2,ts3
	296
	297	! included by RW for temporary comparison
	298	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
	299
	300	! included by RW to compute tracer column densities
	301	real qcol(ngridmx,nqmx)
	302
	303	! included by RW for H2O precipitation
	304	REAL zdtrain(ngridmx,nlayermx)
	305	REAL zdqrain(ngridmx,nlayermx,nqmx)
	306	REAL zdqsrain(ngridmx)
	307	REAL zdqssnow(ngridmx)
	308	REAL rainout(ngridmx)
	309
	310	! included by RW for H2O Manabe scheme
	311	REAL dtmoist(ngridmx,nlayermx)
	312	REAL dqmoist(ngridmx,nlayermx,nqmx)
	313
	314	real qvap(ngridmx,nlayermx)
	315	real dqvaplscale(ngridmx,nlayermx)
	316	real dqcldlscale(ngridmx,nlayermx)
	317	real cloudfrac(ngridmx,nlayermx)
	318	real rneb_man(ngridmx,nlayermx)
	319	real rneb_lsc(ngridmx,nlayermx)
	320
	321	! included by RW to account for surface cooling by evaporation
	322	REAL dtsurfh2olat(ngridmx)
	323
	324	! included by RW to test energy conservation
	325	REAL dEtot, dEtots, masse, vabs, dvabs
	326
	327	! included by RW to test water conservation
	328	real h2otot
	329
	330	! included by RW to allow variations in cp with location
	331	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
	332	REAL rcp3D(ngridmx,nlayermx) ! R / specific heat capacity
	333	real cppNI, rcpNI, nnu ! last one just for Seb version
	334	REAL zpopskNI(ngridmx,nlayermx)
	335
	336	! included by RW to make 1D saves not every timestep
	337	integer countG1D
	338	save countG1D
	339
	340	c=======================================================================
	341
	342	kastprof=.false.
	343
	344	c 1. Initialisation:
	345	c -----------------
	346
	347	c 1.1 Initialisation only at first call
	348	c ---------------------------------------
	349	IF (firstcall) THEN
	350
	351	if(ngrid.eq.1)then
	352	!saveG1D=day_step*10
	353	saveG1D=1
	354	countG1D=1
	355	endif
	356
	357	c variables set to 0
	358	c ~~~~~~~~~~~~~~~~~~
	359	call zerophys(ngrid*nlayer,dtrad)
	360	call zerophys(ngrid,fluxrad)
	361
	362	c initialize tracer names, indexes and properties
	363	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	364	tracerdyn=tracer
	365	IF (tracer) THEN
	366	CALL initracer()
	367	ENDIF ! end tracer
	368
	369
	370	c read startfi (initial parameters)
	371	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	372
	373	CALL phyetat0 ("startfi.nc",0,0,
	374	& nsoilmx,nq,
	375	& day_ini,time_phys,
	376	& tsurf,tsoil,emis,q2,qsurf)
	377
	378	if (pday.ne.day_ini) then
	379	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	380	& "physics and dynamics"
	381	write(,) "dynamics day: ",pday
	382	write(,) "physics day: ",day_ini
	383	stop
	384	endif
	385
	386	write (,) 'In physiq day_ini =', day_ini
	387
	388
	389	c Initialize albedo and orbital calculation
	390	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	391	CALL surfini(ngrid,qsurf,albedo)
	392	! surely for good CO2 ice feedback this must be called every turn?
	393	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	394
	395
	396	if(tlocked)then
	397	print*,'Planet is tidally locked at resonance n=',nres
	398	print*,'Make sure you have the right rotation rate!!!'
	399	endif
	400
	401
	402	c initialize soil
	403	c ~~~~~~~~~~~~~~~
	404	IF (callsoil) THEN
	405	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	406	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	407	ELSE
	408	PRINT*,'WARNING! Thermal conduction in the soil turned off'
	409	DO ig=1,ngrid
	410	capcal(ig)=1.e6
	411	fluxgrd(ig)=0.
	412	ENDDO
	413	ENDIF
	414	icount=1
	415
	416	c Set temperature just above condensation temperature (for Early Mars)
	417	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	418	if(nearco2cond) then
	419	write(,)' WARNING, starting at Tcond+1K !!!'
	420	DO l=1, nlayer
	421	DO ig=1,ngrid
	422	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4
	423	& -pt(ig,l)) / ptimestep
	424	ENDDO
	425	ENDDO
	426	endif
	427
	428	if(meanOLR)then
	429	! to record global radiative balance
	430	call system('rm -f rad_bal.out')
	431	endif
	432
	433	if(water)then
	434	! initialise variables for water cycle
	435
	436	if(ngrid.eq.1)then
	437	qsurf(1,igcm_h2o_ice) = 100.0
	438	DO l=1, nlayer
	439	pq(1,l,igcm_h2o_ice)=0.0!1.0e-2
	440	enddo
	441	endif
	442
	443	endif
	444	call su_watercycle ! even if we don't have a water cycle, we might
	445	! need epsi for the wvp definitions in callcorrk.F
	446
	447	ENDIF ! (end of "if firstcall")
	448
	449	c ---------------------------------------------------
	450	c 1.2 Initializations done at every physical timestep:
	451	c ---------------------------------------------------
	452	c
	453	IF (ngrid.NE.ngridmx) THEN
	454	PRINT*,'STOP in PHYSIQ'
	455	PRINT*,'Probleme de dimensions :'
	456	PRINT*,'ngrid = ',ngrid
	457	PRINT*,'ngridmx = ',ngridmx
	458	STOP
	459	ENDIF
	460
	461	c Initialize various variables
	462	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	463	call zerophys(ngrid*nlayer, pdv)
	464	call zerophys(ngrid*nlayer, pdu)
	465	if ( (.not.nearco2cond).and.(.not.firstcall) ) then
	466	call zerophys(ngrid*nlayer,pdt)
	467	endif ! this was the source of an evil bug...
	468	call zerophys(ngridnlayernq, pdq)
	469	call zerophys(ngrid, pdpsrf)
	470	call zerophys(ngrid, zflubid)
	471	call zerophys(ngrid, zdtsurf)
	472	call zerophys(ngrid*nq, dqsurf)
	473	igout=ngrid/2+1
	474
	475	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	476
	477	c Compute Solar Longitude (Ls) :
	478	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	479	if (season) then
	480	call solarlong(zday,zls)
	481	else
	482	call solarlong(float(day_ini),zls)
	483	end if
	484
	485	c Compute geopotential between layers
	486	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	487
	488	DO l=1,nlayer
	489	DO ig=1,ngrid
	490	zzlay(ig,l)=pphi(ig,l)/g
	491	ENDDO
	492	ENDDO
	493	DO ig=1,ngrid
	494	zzlev(ig,1)=0.
	495	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	496	ENDDO
	497	DO l=2,nlayer
	498	DO ig=1,ngrid
	499	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	500	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	501	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	502	ENDDO
	503	ENDDO
	504
	505	! Potential temperature calculation not the same in physiq and dynamic...
	506
	507	c Compute potential temperature
	508	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	509	if(nonideal)then
	510
	511	DO l=1,nlayer
	512	DO ig=1,ngrid
	513	call calc_cpp3d(cppNI,rcpNI,pt(ig,l),pplay(ig,l))
	514	cpp3D(ig,l) = cppNI
	515	rcp3D(ig,l) = rcpNI
	516	ENDDO
	517	ENDDO
	518
	519	DO l=1,nlayer
	520	DO ig=1,ngrid
	521
	522	! nnu=0.35
	523	! zh(ig,l) = (pt(ig,l)**nnu +
	524	! & nnu(460nnu)(8.314/1000.)*
	525	! & log(pplev(ig,1)/pplay(ig,l)) )**(-nnu) ! Sebastian's version
	526	! zpopskNI(ig,l) = (pplay(ig,l)/pplev(ig,1))**rcp3D(ig,l)
	527
	528	zh(ig,l) = pt(ig,l)*1.0024 - 0.659
	529	zpopskNI(ig,l) = pt(ig,l)/zh(ig,l)
	530	! we're only after zpopskNI here, not zh
	531	! zh is calculated seperately before both vdifc and convadj
	532	ENDDO
	533	ENDDO
	534
	535	endif
	536	! else
	537
	538	DO l=1,nlayer
	539	DO ig=1,ngrid
	540	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	541	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	542	ENDDO
	543	ENDDO
	544
	545	! endif
	546
	547
	548	c-----------------------------------------------------------------------
	549	c 2. Compute radiative tendencies :
	550	c-----------------------------------------------------------------------
	551
	552	IF (callrad) THEN
	553	IF( MOD(icount-1,iradia).EQ.0) THEN
	554
	555	c Local stellar zenith angle
	556	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	557	CALL orbite(zls,dist_sol,declin)
	558
	559	IF (tlocked) THEN
	560	ztim1=SIN(declin)
	561	ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
	562	ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
	563
	564	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	565	s ztim1,ztim2,ztim3,mu0,fract)
	566
	567	ELSEIF (diurnal) THEN
	568	ztim1=SIN(declin)
	569	ztim2=COS(declin)COS(2.pi*(zday-.5))
	570	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	571
	572	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	573	s ztim1,ztim2,ztim3,mu0,fract)
	574
	575	ELSE
	576
	577	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
	578	! WARNING: this function appears not to work in 1D
	579
	580	ENDIF
	581
	582	c Call main radiative transfer scheme
	583	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	584
	585	c Radiative transfer
	586	c ------------------
	587
	588	if (corrk) then
	589	if(kastprof)then
	590	call kastprof_fn(tsurf,pt,pplay,pplev)
	591	print*,'pt',pt
	592	print*,'pplay',pplay
	593	endif
	594
	595	CALL callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
	596	& albedo,emis,mu0,pplev,pplay,pt,
	597	& tsurf,fract,dist_sol,igout,aerosol,cpp3D,
	598	& zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
	599	& fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,
	600	& firstcall,lastcall)
	601
	602
	603
	604	c Radiative flux from the sky absorbed by the surface (W.m-2)
	605	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	606	DO ig=1,ngrid
	607	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	608	& +fluxsurf_sw(ig)*(1.-albedo(ig))
	609	ENDDO
	610
	611	c Net atmospheric radiative heating rate (K.s-1)
	612	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	613
	614	DO l=1,nlayer
	615	DO ig=1,ngrid
	616	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	617	ENDDO
	618	ENDDO
	619
	620	else
	621
	622	c Atmosphere has no radiative effect
	623	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	624	DO ig=1,ngrid
	625	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_sol**2
	626	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
	627	fluxtop_sw(ig) = fluxtop_dn(ig)*albedo(ig)
	628	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
	629	ENDDO ! radiation skips the atmosphere entirely
	630
	631
	632	DO l=1,nlayer
	633	DO ig=1,ngrid
	634	dtrad(ig,l)=0.0
	635	ENDDO
	636	ENDDO ! hence no atmospheric radiative heating
	637
	638	endif ! if corrk
	639
	640	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	641
	642	! Transformation of the radiative tendencies:
	643	! -------------------------------------------
	644
	645	! Net radiative surface flux (W.m-2)
	646	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	647
	648	DO ig=1,ngrid
	649	zplanck(ig)=tsurf(ig)*tsurf(ig)
	650	zplanck(ig)=emis(ig)*
	651	& stephanzplanck(ig)zplanck(ig)
	652	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
	653	ENDDO
	654
	655	DO l=1,nlayer
	656	DO ig=1,ngrid
	657	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	658	ENDDO
	659	ENDDO
	660
	661	!-------------------------
	662	! test energy conservation
	663	if(enertest)then
	664	dEtot=0.0
	665	dEtots=0.0
	666	DO ig = 1, ngrid
	667	DO l = 1, nlayer
	668	masse = (pplev(ig,l) - pplev(ig,l+1))/g
	669	dEtot = dEtot + cppmassedtrad(ig,l)*area(ig)
	670	!print,'l=',l,'dEtot=',cppmassedtrad(ig,l)area(ig)
	671	ENDDO
	672	dEtots = dEtots + fluxrad(ig)*area(ig)
	673	ENDDO
	674	dEtot=dEtot/totarea
	675	dEtots=dEtots/totarea
	676	print*,'In corrk atmospheric energy change=',dEtot,' W m-2'
	677	print*,'In corrk surface energy change=',dEtots,' W m-2'
	678	endif
	679	!-------------------------
	680
	681	ENDIF ! of IF (callrad)
	682
	683
	684	!-----------------------------------------------------------------------
	685	! 4. Vertical diffusion (turbulent mixing):
	686	! -----------------------------------------
	687
	688	IF (calldifv) THEN
	689
	690	DO ig=1,ngrid
	691	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	692	ENDDO
	693
	694	CALL zerophys(ngrid*nlayer,zdum1)
	695	CALL zerophys(ngrid*nlayer,zdum2)
	696	do l=1,nlayer
	697	do ig=1,ngrid
	698	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	699	enddo
	700	enddo
	701
	702	c Calling vdif (Martian version WITH CO2 condensation)
	703	CALL vdifc(ngrid,nlayer,nq,zpopsk,
	704	& ptimestep,capcal,lwrite,
	705	& pplay,pplev,zzlay,zzlev,z0,
	706	& pu,pv,zh,pq,tsurf,emis,qsurf,
	707	& zdum1,zdum2,zdh,pdq,zflubid,
	708	& zdudif,zdvdif,zdhdif,zdtsdif,q2,
	709	& zdqdif,zdqsdif)
	710
	711	DO l=1,nlayer
	712	DO ig=1,ngrid
	713	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	714	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	715	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	716
	717	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	718	ENDDO
	719	ENDDO
	720
	721	DO ig=1,ngrid
	722	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	723	ENDDO
	724
	725	if (tracer) then
	726	DO iq=1, nq
	727	DO l=1,nlayer
	728	DO ig=1,ngrid
	729	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	730	ENDDO
	731	ENDDO
	732	ENDDO
	733	DO iq=1, nq
	734	DO ig=1,ngrid
	735	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	736	ENDDO
	737	ENDDO
	738
	739	end if ! of if (tracer)
	740
	741	!print*,'surface temp loss:'
	742	!print*, dtsurfh2olat
	743
	744	!-------------------------
	745	! test energy conservation
	746	if(enertest)then
	747	dEtot=0.0
	748	dEtots=0.0
	749	DO ig = 1, ngrid
	750	DO l = 1, nlayer
	751	masse = (pplev(ig,l) - pplev(ig,l+1))/g
	752
	753	dEtot = dEtot + cppmassezdtdif(ig,l)*area(ig)
	754
	755	vabs = sqrt(pdu(ig,l)2 + pdv(ig,l)2)
	756	dvabs = sqrt(zdudif(ig,l)2 + zdvdif(ig,l)2)
	757	dEtot = dEtot + massevabsdvabs*area(ig)
	758	ENDDO
	759	dEtots = dEtots + capcal(ig)zdtsdif(ig)area(ig)
	760	ENDDO
	761	dEtot=dEtot/totarea
	762	dEtots=dEtots/totarea
	763	print*,'In difv atmospheric energy change=',dEtot,' W m-2'
	764	print*,'In difv surface energy change=',dEtots,' W m-2'
	765	endif
	766	!-------------------------
	767
	768	ELSE
	769
	770	DO ig=1,ngrid
	771	zdtsurf(ig)=zdtsurf(ig)+
	772	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	773	ENDDO
	774
	775	ENDIF ! of IF (calldifv)
	776
	777
	778	!-----------------------------------------------------------------------
	779	! 5. Dry convective adjustment:
	780	! -----------------------------
	781
	782	IF(calladj) THEN
	783
	784	DO l=1,nlayer
	785	DO ig=1,ngrid
	786	if(nonideal)then
	787	zdh(ig,l)=pdt(ig,l)/zpopskNI(ig,l)
	788	else
	789	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	790	endif
	791	ENDDO
	792	ENDDO
	793	CALL zerophys(ngrid*nlayer,zduadj)
	794	CALL zerophys(ngrid*nlayer,zdvadj)
	795	CALL zerophys(ngrid*nlayer,zdhadj)
	796
	797	if(nonideal)then
	798	CALL convadj(ngrid,nlayer,nq,ptimestep,
	799	& pplay,pplev,zpopskNI,
	800	& pu,pv,zh,pq,
	801	& pdu,pdv,zdh,pdq,
	802	& zduadj,zdvadj,zdhadj,
	803	& zdqadj)
	804	else
	805	CALL convadj(ngrid,nlayer,nq,ptimestep,
	806	& pplay,pplev,zpopsk,
	807	& pu,pv,zh,pq,
	808	& pdu,pdv,zdh,pdq,
	809	& zduadj,zdvadj,zdhadj,
	810	& zdqadj)
	811	endif
	812
	813
	814	DO l=1,nlayer
	815	DO ig=1,ngrid
	816	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	817	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	818	if(nonideal)then
	819	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopskNI(ig,l)
	820	zdtadj(ig,l)=zdhadj(ig,l)*zpopskNI(ig,l) ! for diagnostic only
	821	else
	822	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	823	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	824	endif
	825	ENDDO
	826	ENDDO
	827
	828	if(tracer) then
	829	DO iq=1, nq
	830	DO l=1,nlayer
	831	DO ig=1,ngrid
	832	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	833	ENDDO
	834	ENDDO
	835	ENDDO
	836	end if
	837
	838	!-------------------------
	839	! test energy conservation
	840	if(enertest)then
	841	dEtot=0.0
	842	DO ig = 1, ngrid
	843	DO l = 1, nlayer
	844	masse = (pplev(ig,l) - pplev(ig,l+1))/g
	845	dEtot = dEtot + cppmassezdtadj(ig,l)*area(ig)
	846	ENDDO
	847	ENDDO
	848	dEtot=dEtot/totarea
	849	print*,'In convadj atmospheric energy change=',dEtot,' W m-2'
	850	endif
	851	!-------------------------
	852
	853	ENDIF ! of IF(calladj)
	854
	855	!-----------------------------------------------------------------------
	856	! 6. Carbon dioxide condensation-sublimation:
	857	! -------------------------------------------
	858
	859	IF (co2cond) THEN
	860	if (.not.tracer) then
	861	print*,'We need a CO2 ice tracer to condense CO2'
	862	call abort
	863	endif
	864
	865	call condens_co2cloud(ngrid,nlayer,nq,ptimestep,
	866	& capcal,pplay,pplev,tsurf,pt,
	867	& pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	868	& qsurf(1,igcm_co2_ice),albedo,emis,
	869	& zdtc,zdtsurfc,pdpsrf,zduc,zdvc,
	870	& zdqc,reffrad,cpp3D)
	871
	872	DO l=1,nlayer
	873	DO ig=1,ngrid
	874	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	875	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	876	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	877	ENDDO
	878	ENDDO
	879	DO ig=1,ngrid
	880	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	881	ENDDO
	882
	883	! IF (tracer) THEN
	884	DO iq=1, nq ! should use new notation here !
	885	DO l=1,nlayer
	886	DO ig=1,ngrid
	887	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	888	ENDDO
	889	ENDDO
	890	ENDDO
	891	!NB: we do not add tendency on surface co2ice ;
	892	! since qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
	893	! ENDIF ! of IF (tracer)
	894
	895
	896	!-------------------------
	897	! test energy conservation
	898	if(enertest)then
	899	dEtot=0.0
	900	DO ig = 1, ngrid
	901	DO l = 1, nlayer
	902	masse = (pplev(ig,l) - pplev(ig,l+1))/g
	903	dEtot = dEtot + cppmassezdtc(ig,l)
	904	ENDDO
	905	dEtot = dEtot + capcal(ig)*zdtsurfc(ig)
	906	ENDDO
	907	print*,'In co2cloud energy change=',dEtot,' W m-2'
	908	endif ! we're missing the latent heat...
	909	!-------------------------
	910
	911	ENDIF ! of IF (co2cond)
	912
	913
	914
	915	c-----------------------------------------------------------------------
	916	c 7. Specific parameterizations for tracers
	917	c -----------------------------------------
	918
	919	if (tracer) then
	920
	921	c 7a. Water and ice
	922	c ---------------
	923
	924	c ---------------------------------------
	925	c Water ice condensation in the atmosphere
	926	c ----------------------------------------
	927	IF (water) THEN
	928
	929	if(watercond)then
	930
	931	if(1.eq.2)then
	932	call moistadj(pt,pq,pplev,pplay,dtmoist,dqmoist,
	933	& ptimestep,rneb_man)
	934
	935	DO l=1,nlayer
	936	DO ig=1,ngrid
	937	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	938	& dqmoist(ig,l,igcm_h2o_vap)
	939	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	940	& dqmoist(ig,l,igcm_h2o_ice)
	941	pdt(ig,l)=pdt(ig,l)+dtmoist(ig,l)
	942	ENDDO
	943	ENDDO
	944	endif ! moistadj turned off for now
	945
	946	DO l=1,nlayer
	947	DO ig=1,ngrid
	948	qvap(ig,l)=pq(ig,l,igcm_h2o_vap)
	949	ENDDO
	950	ENDDO
	951	call largescale2(ptimestep,pplev,pplay,pt,qvap,
	952	s zdtlscale, dqvaplscale,dqcldlscale,rneb_lsc)
	953	DO l=1,nlayer
	954	DO ig=1,ngrid
	955	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	956	& dqvaplscale(ig,l)
	957	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	958	& dqcldlscale(ig,l)
	959	! pdt(ig,l)=pdt(ig,l)+zdtlscale(ig,l)
	960	! for the moment this is removed, as it dont conserve energy
	961	ENDDO
	962	ENDDO
	963
	964	! compute cloud fraction
	965	DO l = 1, nlayer
	966	DO ig = 1,ngrid
	967	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
	968	ENDDO
	969	ENDDO
	970
	971	endif ! of IF (watercondense)
	972
	973	if(waterrain)then
	974
	975	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,
	976	& zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
	977
	978	DO l=1,nlayer
	979	DO ig=1,ngrid
	980	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	981	& zdqrain(ig,l,igcm_h2o_vap)
	982	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	983	& zdqrain(ig,l,igcm_h2o_ice)
	984	!pdt(ig,l)=pdt(ig,l)+zdtrain(ig,l)
	985	ENDDO
	986	ENDDO
	987
	988	DO ig=1,ngrid
	989	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
	990	& zdqsrain(ig)+zdqssnow(ig)
	991	rainout(ig)=zdqsrain(ig)+zdqssnow(ig) ! diagnostic
	992	!rainout(ig)=zdqssnow(ig)
	993	ENDDO
	994
	995	end if ! of IF (waterrain)
	996
	997	end if ! of IF (water)
	998
	999
	1000	c 7c. Aerosol particles
	1001	c -------------------
	1002
	1003	c -------------
	1004	c Sedimentation
	1005	c -------------
	1006	IF (sedimentation) THEN
	1007	call zerophys(ngridnlayernq, zdqsed)
	1008	call zerophys(ngrid*nq, zdqssed)
	1009
	1010	call callsedim(ngrid,nlayer, ptimestep,
	1011	& pplev,zzlev, pt,
	1012	& pq, pdq, zdqsed, zdqssed,nq)
	1013
	1014	DO iq=2, nq ! should be updated
	1015	DO l=1,nlayer
	1016	DO ig=1,ngrid
	1017	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1018	ENDDO
	1019	ENDDO
	1020	ENDDO
	1021	DO iq=2, nq
	1022	DO ig=1,ngrid
	1023	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1024	ENDDO
	1025	ENDDO
	1026	END IF ! of IF (sedimentation)
	1027
	1028
	1029	c 7d. Updates
	1030	c ---------
	1031
	1032	DO iq=1, nq
	1033	DO ig=1,ngrid
	1034
	1035	c ---------------------------------
	1036	c Updating tracer budget on surface
	1037	c ---------------------------------
	1038	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1039
	1040	ENDDO ! (ig)
	1041	ENDDO ! (iq)
	1042
	1043	endif ! of if (tracer)
	1044
	1045
	1046	!-----------------------------------------------------------------------
	1047	! 9. Surface and sub-surface soil temperature
	1048	!-----------------------------------------------------------------------
	1049	!
	1050	!
	1051	! 9.1 Increment Surface temperature:
	1052	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1053
	1054	DO ig=1,ngrid
	1055	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1056	ENDDO
	1057
	1058	!
	1059	! 9.2 Compute soil temperatures and subsurface heat flux:
	1060	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1061	IF (callsoil) THEN
	1062	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
	1063	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1064	ENDIF
	1065
	1066	!-------------------------
	1067	! test energy conservation
	1068	if(enertest)then
	1069	dEtots=0.0
	1070	DO ig = 1, ngrid
	1071	dEtots = dEtots + capcal(ig)zdtsurf(ig)area(ig)
	1072	ENDDO
	1073	dEtots=dEtots/totarea
	1074	print*,'Surface energy change=',dEtots,' W m-2'
	1075	endif
	1076	!-------------------------
	1077
	1078
	1079
	1080
	1081	!-----------------------------------------------------------------------
	1082	! 10. Perform diagnostics and write output files
	1083	!-----------------------------------------------------------------------
	1084
	1085	! ---------------------------------------------------------
	1086	! Check the energy balance of the simulation during the run
	1087	! ---------------------------------------------------------
	1088	flux1 = 0
	1089	flux2 = 0
	1090	flux3 = 0
	1091	do ig=1,ngrid
	1092
	1093	flux1 = flux1 +
	1094	& area(ig)*(fluxtop_dn(ig) - fluxtop_sw(ig))
	1095	flux2 = flux2 + area(ig)*fluxtop_lw(ig)
	1096	flux3 = flux3 + area(ig)*fluxtop_dn(ig)
	1097	fluxabs_sw(ig)=fluxtop_dn(ig) - fluxtop_sw(ig)
	1098
	1099	end do
	1100	print*,'Incident solar flux, absorbed solar flux, OLR (W/m2)'
	1101	print*, flux3/totarea, ' ', flux1/totarea ,
	1102	& ' = ', flux2/totarea
	1103
	1104	if(meanOLR)then
	1105	! to record global radiative balance
	1106	open(92,file="rad_bal.out",form='formatted',access='append')
	1107	write(92,*) zday,flux3/totarea,flux1/totarea,flux2/totarea
	1108	close(92)
	1109	endif
	1110
	1111	if(kastprof)then ! for f(p,T) computation (not needed in universal model)
	1112	open(15,file='ene_bal.out',form='formatted',access='append')
	1113	write(15,*) flux1
	1114	close(15)
	1115	stop
	1116	endif
	1117
	1118	! Surface temperature information
	1119	ts1 = 0
	1120	ts2 = 999
	1121	ts3 = 0
	1122	do ig=1,ngrid
	1123	ts1 = ts1 + area(ig)*tsurf(ig)
	1124	ts2 = min(ts2,tsurf(ig))
	1125	ts3 = max(ts3,tsurf(ig))
	1126	end do
	1127	print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2,
	1128	& ' Max Tsurf =',ts3
	1129
	1130	! -------------------------------
	1131	! Dynamical fields incrementation
	1132	! -------------------------------
	1133	! (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1134	! temperature, zonal and meridional wind
	1135	DO l=1,nlayer
	1136	DO ig=1,ngrid
	1137	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1138	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1139	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1140
	1141	! diagnostic
	1142	zdtdyn(ig,l)=ztprevious(ig,l)-pt(ig,l)
	1143	ztprevious(ig,l)=zt(ig,l)
	1144	ENDDO
	1145	ENDDO
	1146	if(firstcall)call zerophys(zdtdyn)
	1147
	1148	! dynamical heating diagnostic
	1149	DO ig=1,ngrid
	1150	fluxdyn(ig)=0.
	1151	DO l=1,nlayer
	1152	if(nonideal)then
	1153	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
	1154	& (pplev(ig,l)-pplev(ig,l+1))cpp3D(ig,l)/g
	1155	else
	1156	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
	1157	& (pplev(ig,l)-pplev(ig,l+1))cpp/g
	1158	endif
	1159	ENDDO
	1160	ENDDO
	1161
	1162
	1163	! tracers
	1164	DO iq=1, nq
	1165	DO l=1,nlayer
	1166	DO ig=1,ngrid
	1167	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1168	ENDDO
	1169	ENDDO
	1170	ENDDO
	1171
	1172	! surface pressure
	1173	DO ig=1,ngrid
	1174	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1175	ENDDO
	1176
	1177	! pressure
	1178	DO l=1,nlayer
	1179	DO ig=1,ngrid
	1180	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1181	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1182	ENDDO
	1183	ENDDO
	1184
	1185
	1186	! ---------------------------------------------------------
	1187	! Compute column amounts (kg m-2) if tracers are enabled
	1188	! ---------------------------------------------------------
	1189	if(tracer)then
	1190	call zerophys(ngrid*nq,qcol)
	1191	do iq=1,nq
	1192	DO ig=1,ngrid
	1193	DO l=1,nlayer
	1194	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *
	1195	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1196	enddo
	1197	enddo
	1198	enddo
	1199	endif
	1200
	1201
	1202	! ---------------------------------------------------------
	1203	! Test for water conservation if water is enabled
	1204	! ---------------------------------------------------------
	1205	if(water)then
	1206
	1207	h2otot = 0
	1208	do ig=1,ngrid
	1209	h2otot = h2otot + area(ig)*
	1210	& (qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap)
	1211	& +qsurf(ig,igcm_h2o_ice)+qsurf(ig,igcm_h2o_vap))
	1212	end do
	1213	print*,'Total water amount (kg): ',h2otot!/totarea
	1214	endif
	1215
	1216
	1217	! TEMPORARY: save water cycle diagnostics
	1218	! IF((ngrid.eq.1).and.lastcall) THEN
	1219	! open(13,file='Ts.out')
	1220	! open(14,file='T.out')
	1221	! open(15,file='p.out')
	1222	! open(16,file='vap.out')
	1223	! open(17,file='alt.out')
	1224	! ig=1
	1225	! write(13,*) tsurf(ig)
	1226	! DO l=1,nlayer
	1227	! write(14,*) pt(ig,l)
	1228	! write(15,*) pplay(ig,l)
	1229	! write(16,*) pq(ig,l,igcm_h2o_vap)
	1230	! write(17,*) zzlay(ig,l)
	1231	! enddo
	1232	! close(13)
	1233	! close(14)
	1234	! close(15)
	1235	! close(16)
	1236	! close(17)
	1237	! ENDIF
	1238
	1239
	1240	! TEMPORARY: save data for LBL
	1241	! IF((ngrid.eq.1).and.lastcall) THEN
	1242	! open(14,file='T.dat')
	1243	! open(15,file='p.dat')
	1244	! ig=1
	1245	! write(14,*) nlayer+1
	1246	! write(15,*) nlayer+1
	1247	! write(14,*) tsurf(ig)
	1248	! write(15,*) pplev(ig,1)
	1249	! DO l=1,nlayer
	1250	! write(14,*) pt(ig,l)
	1251	! write(15,*) pplay(ig,l)
	1252	! enddo
	1253	! close(14)
	1254	! close(15)
	1255	! ENDIF
	1256
	1257
	1258	IF (ngrid.NE.1) THEN
	1259	print,'Ls =',zls180./pi
	1260
	1261	! -------------------------------------------------------------------
	1262	! Writing NetCDF file "RESTARTFI" at the end of the run
	1263	! -------------------------------------------------------------------
	1264	! Note: 'restartfi' is stored just before dynamics are stored
	1265	! in 'restart'. Between now and the writting of 'restart',
	1266	! there will have been the itau=itau+1 instruction and
	1267	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1268	! thus we store for time=time+dtvr
	1269
	1270	IF(lastcall) THEN
	1271	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1272
	1273	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1274	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1275	. ptimestep,pday,
	1276	. ztime_fin,tsurf,tsoil,emis,q2,qsurf,
	1277	. area,albedodat,inertiedat,zmea,zstd,zsig,
	1278	. zgam,zthe)
	1279	ENDIF
	1280
	1281	! -----------------------------------------------------------------
	1282	! Saving statistics :
	1283	! -----------------------------------------------------------------
	1284	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1285	! which can later be used to make the statistic files of the run:
	1286	! "stats") only possible in 3D runs !
	1287
	1288
	1289	IF (callstats) THEN
	1290
	1291	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1292	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1293	call wstats(ngrid,"fluxsurf_lw",
	1294	. "Thermal IR radiative flux to surface","W.m-2",2,
	1295	. fluxsurf_lw)
	1296	call wstats(ngrid,"fluxsurf_sw",
	1297	. "Solar radiative flux to surface","W.m-2",2,
	1298	. fluxsurf_sw_tot)
	1299	call wstats(ngrid,"fluxtop_lw",
	1300	. "Thermal IR radiative flux to space","W.m-2",2,
	1301	. fluxtop_lw)
	1302	call wstats(ngrid,"fluxtop_sw",
	1303	. "Solar radiative flux to space","W.m-2",2,
	1304	. fluxtop_sw_tot)
	1305
	1306	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1307	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1308	call wstats(ngrid,"v","Meridional (North-South) wind",
	1309	. "m.s-1",3,zv)
	1310	call wstats(ngrid,"w","Vertical (down-up) wind",
	1311	. "m.s-1",3,pw)
	1312	call wstats(ngrid,"q2",
	1313	. "Boundary layer eddy kinetic energy","m2.s-2",3,q2)
	1314
	1315	if (tracer) then
	1316	if (water) then
	1317	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1318	& *mugaz/mmol(igcm_h2o_vap)
	1319	call wstats(ngrid,"vmr_h2ovapor",
	1320	. "H2O vapor volume mixing ratio","mol/mol",
	1321	. 3,vmr)
	1322	endif ! of if (water)
	1323
	1324	endif !tracer
	1325
	1326	IF(lastcall) THEN
	1327	write (,) "Writing stats..."
	1328	call mkstats(ierr)
	1329	ENDIF
	1330	ENDIF !if callstats
	1331
	1332
	1333
	1334	! ----------------------------------------------------------------------
	1335	! output in netcdf file "DIAGFI", containing any variable for diagnostic
	1336	! (output with period "ecritphy", set in "run.def")
	1337	! ----------------------------------------------------------------------
	1338	! WRITEDIAGFI can ALSO be called from any other subroutine
	1339	! for any variable!
	1340
	1341	! is it not preferable to keep all the calls in one place?
	1342
	1343	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1344	& tsurf)
	1345	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1346	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1347	! call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
	1348	! call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1349	! call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1350	! call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1351	! call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
	1352	! call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
	1353	! call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh)
	1354	! call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay)
	1355
	1356	! Total energy balance diagnostics
	1357	call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
	1358	& 2,fluxtop_dn)
	1359	call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
	1360	& 2,fluxabs_sw)
	1361	call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
	1362	& 2,fluxtop_lw)
	1363	call WRITEDIAGFI(ngrid,"DYN","dynamical heat input","W m-2",
	1364	& 2,fluxdyn)
	1365
	1366	! Temporary inclusions for heating diagnostics
	1367	! call WRITEDIAGFI(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
	1368	! call WRITEDIAGFI(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	1369	! call WRITEDIAGFI(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	1370
	1371	! Output tracers
	1372	if (tracer) then
	1373	do iq=1,nq
	1374	! call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	1375	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
	1376	& trim(noms(iq))//'_surf',
	1377	& 'kg/kg',2,qsurf(1,iq) )
	1378
	1379	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
	1380	& trim(noms(iq))//'_col',
	1381	& 'kg m^-2',2,qcol(1,iq) )
	1382
	1383
	1384
	1385	if(waterrain)then
	1386	CALL WRITEDIAGFI(ngridmx,"rain",
	1387	& "total precipitation","kg m-2 s-1",2,rainout)
	1388	endif
	1389
	1390	! call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_aero',
	1391	! & trim(noms(iq))//'_aero',
	1392	! & 'kg/kg',3,aerosol(1,1,iq))
	1393
	1394	call WRITEDIAGFI(ngridmx,"tau_col",
	1395	& "Total aerosol optical depth","[]",2,tau_col)
	1396
	1397	! call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_reff',
	1398	! & trim(noms(iq))//'_reff',
	1399	! & 'm',3,reffrad(1,1,iq))
	1400
	1401	enddo
	1402	endif
	1403
	1404	ELSE ! if(ngrid.eq.1)
	1405
	1406
	1407	! ----------------------------------------------------------------------
	1408	! Output in grads file "g1d" (ONLY when using testphys1d)
	1409	! ----------------------------------------------------------------------
	1410
	1411	if(countG1D.eq.saveG1D)then
	1412
	1413	! call WRITEDIAGFI(ngrid,"tsurf","Surface temperature",
	1414	! & "K",0,tsurf)
	1415	CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	1416	CALL writeg1d(ngrid,1,ps,'ps','Pa')
	1417	CALL writeg1d(ngrid,nlayer,zt,'T','K')
	1418	CALL writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg')
	1419	CALL writeg1d(ngrid,nlayer,aerosol,'aerosol','SI')
	1420	CALL writeg1d(ngrid,nlayer,reffrad,'reffrad','SI')
	1421	CALL writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI')
	1422	CALL writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI')
	1423	CALL writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI')
	1424
	1425	! radiation balance (optional)
	1426	CALL writeg1d(ngrid,1,flux3/totarea,'ISR','W m-2')
	1427	CALL writeg1d(ngrid,1,flux1/totarea,'ASR','W m-2')
	1428	CALL writeg1d(ngrid,1,flux2/totarea,'OLR','W m-2')
	1429
	1430	if(tracer) then
	1431	do iq=1,nq
	1432	CALL writeg1d(ngrid,1,qsurf(1,iq),
	1433	& trim(noms(iq))//'_s','kg m^-2')
	1434	CALL writeg1d(ngrid,1,qcol(1,iq),
	1435	& trim(noms(iq))//'_c','kg m^-2')
	1436	CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	1437	end do
	1438
	1439	if(waterrain)then
	1440	CALL writeg1d(ngrid,1,rainout,
	1441	& 'rainfall','kg m-2 s-1')
	1442	endif
	1443	end if
	1444
	1445	countG1D=1
	1446	else
	1447	countG1D=countG1D+1
	1448	endif ! if time to save
	1449
	1450	END IF ! if(ngrid.ne.1)
	1451
	1452	icount=icount+1
	1453
	1454	RETURN
	1455	END

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

Download in other formats:

Original Format