Context Navigation

physiq.F90 @ 586

Last change on this file since 586 was 586, checked in by jleconte, 13 years ago
removed cpp3D and nonideal stuff.
Property svn:executable set to ``*
File size: 75.1 KB

Line
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,L_NSPECTI,L_NSPECTV
10	use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
11	use gases_h
12	implicit none
13
14
15	!==================================================================
16	!
17	! Purpose
18	! -------
19	! Central subroutine for all the physics parameterisations in the
20	! universal model. Originally adapted from the Mars LMDZ model.
21	!
22	! The model can be run without or with tracer transport
23	! depending on the value of "tracer" in file "callphys.def".
24	!
25	!
26	! It includes:
27	!
28	! 1. Initialization:
29	! 1.1 Firstcall initializations
30	! 1.2 Initialization for every call to physiq
31	! 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
32	! 2. Compute radiative transfer tendencies
33	! (longwave and shortwave).
34	! 4. Vertical diffusion (turbulent mixing):
35	! 5. Convective adjustment
36	! 6. Condensation and sublimation of gases (currently just CO2).
37	! 7. TRACERS
38	! 7a. water and water ice
39	! 7c. other schemes for tracer transport (lifting, sedimentation)
40	! 7d. updates (pressure variations, surface budget)
41	! 9. Surface and sub-surface temperature calculations
42	! 10. Write outputs :
43	! - "startfi", "histfi" if it's time
44	! - Saving statistics if "callstats = .true."
45	! - Output any needed variables in "diagfi"
46	! 10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
47	!
48	! arguments
49	! ---------
50	!
51	! input
52	! -----
53	! ecri period (in dynamical timestep) to write output
54	! ngrid Size of the horizontal grid.
55	! All internal loops are performed on that grid.
56	! nlayer Number of vertical layers.
57	! nq Number of advected fields
58	! firstcall True at the first call
59	! lastcall True at the last call
60	! pday Number of days counted from the North. Spring
61	! equinoxe.
62	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
63	! ptimestep timestep (s)
64	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
65	! pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
66	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
67	! pu(ngrid,nlayer) u component of the wind (ms-1)
68	! pv(ngrid,nlayer) v component of the wind (ms-1)
69	! pt(ngrid,nlayer) Temperature (K)
70	! pq(ngrid,nlayer,nq) Advected fields
71	! pudyn(ngrid,nlayer) \
72	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
73	! ptdyn(ngrid,nlayer) / corresponding variables
74	! pqdyn(ngrid,nlayer,nq) /
75	! pw(ngrid,?) vertical velocity
76	!
77	! output
78	! ------
79	!
80	! pdu(ngrid,nlayermx) \
81	! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
82	! pdt(ngrid,nlayermx) / variables due to physical processes.
83	! pdq(ngrid,nlayermx) /
84	! pdpsrf(ngrid) /
85	! tracerdyn call tracer in dynamical part of GCM ?
86	!
87	!
88	! Authors
89	! -------
90	! Frederic Hourdin 15/10/93
91	! Francois Forget 1994
92	! Christophe Hourdin 02/1997
93	! Subroutine completely rewritten by F. Forget (01/2000)
94	! Water ice clouds: Franck Montmessin (update 06/2003)
95	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
96	! New correlated-k radiative scheme: R. Wordsworth (2009)
97	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
98	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
99	! To F90: R. Wordsworth (2010)
100	!
101	!==================================================================
102
103
104	! 0. Declarations :
105	! ------------------
106
107	#include "dimensions.h"
108	#include "dimphys.h"
109	#include "comgeomfi.h"
110	#include "surfdat.h"
111	#include "comsoil.h"
112	#include "comdiurn.h"
113	#include "callkeys.h"
114	#include "comcstfi.h"
115	#include "planete.h"
116	#include "comsaison.h"
117	#include "control.h"
118	#include "tracer.h"
119	#include "watercap.h"
120	#include "netcdf.inc"
121
122	! Arguments :
123	! -----------
124
125	! inputs:
126	! -------
127	integer ngrid,nlayer,nq
128	real ptimestep
129	real pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
130	real pphi(ngridmx,nlayer)
131	real pu(ngridmx,nlayer),pv(ngridmx,nlayer)
132	real pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
133	real pw(ngridmx,nlayer) ! pvervel transmitted by dyn3d
134	real zh(ngridmx,nlayermx) ! potential temperature (K)
135	logical firstcall,lastcall
136
137	real pday
138	real ptime
139	logical tracerdyn
140
141	! outputs:
142	! --------
143	! physical tendencies
144	real pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
145	real pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
146	real pdpsrf(ngridmx) ! surface pressure tendency
147
148
149	! Local saved variables:
150	! ----------------------
151	! aerosol (dust or ice) extinction optical depth at reference wavelength
152	! "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
153	! aerosol optical properties:
154	! real aerosol(ngridmx,nlayermx,naerkind)
155	! this is now internal to callcorrk and hence no longer needed here
156
157	integer day_ini ! Initial date of the run (sol since Ls=0)
158	integer icount ! counter of calls to physiq during the run.
159	real tsurf(ngridmx) ! Surface temperature (K)
160	real tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
161	real albedo(ngridmx) ! Surface albedo
162
163	real albedo0(ngridmx) ! Surface albedo
164	integer rnat(ngridmx) ! added by BC
165	save rnat
166
167	real emis(ngridmx) ! Thermal IR surface emissivity
168	real dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
169	real fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
170	real fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
171	real capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
172	real fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
173	real qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
174	real q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
175
176	save day_ini, icount
177	save tsurf,tsoil
178	save albedo0,albedo,emis,q2
179	save capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
180
181	real, parameter :: stephan = 5.67e-08 ! Stephan Boltzman constant (should really be elsewhere)
182
183	! Local variables :
184	! -----------------
185
186	! aerosol (dust or ice) extinction optical depth at reference wavelength
187	! for the "naerkind" optically active aerosols:
188	real aerosol(ngridmx,nlayermx,naerkind)
189
190	character*80 fichier
191	integer l,ig,ierr,iq,i, tapphys,nw
192
193	real fluxsurf_lw(ngridmx) ! incident LW (IR) surface flux (W.m-2)
194	real fluxsurf_sw(ngridmx) ! incident SW (stellar) surface flux (W.m-2)
195	real fluxtop_lw(ngridmx) ! Outgoing LW (IR) flux to space (W.m-2)
196	real fluxabs_sw(ngridmx) ! Absorbed SW (stellar) flux (W.m-2)
197
198	real fluxtop_dn(ngridmx)
199	real fluxdyn(ngridmx) ! horizontal heat transport by dynamics
200	real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
201	real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
202	save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
203
204
205	real zls ! solar longitude (rad)
206	real zday ! date (time since Ls=0, in martian days)
207	real zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
208	real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
209	real latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
210
211	real reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
212
213	! Tendencies due to various processes:
214	real dqsurf(ngridmx,nqmx)
215	real zdtlw(ngridmx,nlayermx) ! (K/s)
216	real zdtsw(ngridmx,nlayermx) ! (K/s)
217	save zdtlw, zdtsw
218
219	real cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear areas
220	real cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear areas
221	real zdtsurf(ngridmx) ! (K/s)
222	real dtlscale(ngridmx,nlayermx)
223	real zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
224	real zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
225	real zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
226	real zdhadj(ngridmx,nlayermx) ! (K/s)
227	real zdtgw(ngridmx,nlayermx) ! (K/s)
228	real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
229	real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
230	real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
231
232	real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
233	real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
234	real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
235	real zdqadj(ngridmx,nlayermx,nqmx)
236	real zdqc(ngridmx,nlayermx,nqmx)
237	real zdqlscale(ngridmx,nlayermx,nqmx)
238	real zdqslscale(ngridmx,nqmx)
239	real zdqchim(ngridmx,nlayermx,nqmx)
240	real zdqschim(ngridmx,nqmx)
241
242	real zdteuv(ngridmx,nlayermx) ! (K/s)
243	real zdtconduc(ngridmx,nlayermx) ! (K/s)
244	real zdumolvis(ngridmx,nlayermx)
245	real zdvmolvis(ngridmx,nlayermx)
246	real zdqmoldiff(ngridmx,nlayermx,nqmx)
247
248	! Local variables for local calculations:
249	real zflubid(ngridmx)
250	real zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
251	real zdum1(ngridmx,nlayermx)
252	real zdum2(ngridmx,nlayermx)
253	real ztim1,ztim2,ztim3, z1,z2
254	real ztime_fin
255	real zdh(ngridmx,nlayermx)
256	integer length
257	parameter (length=100)
258
259	! local variables only used for diagnostics (output in file "diagfi" or "stats")
260	! ------------------------------------------------------------------------------
261	real ps(ngridmx), zt(ngridmx,nlayermx)
262	real zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
263	real zq(ngridmx,nlayermx,nqmx)
264	character*2 str2
265	character*5 str5
266	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
267	real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
268	save ztprevious
269	real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
270	real qtot1,qtot2 ! total aerosol mass
271	integer igmin, lmin
272	logical tdiag
273
274	real zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
275	real zstress(ngrid), cd
276	real hco2(nqmx), tmean, zlocal(nlayermx)
277	real vmr(ngridmx,nlayermx) ! volume mixing ratio
278
279	real time_phys
280
281	! reinstated by RW for diagnostic
282	real tau_col(ngridmx)
283	save tau_col
284
285	! included by RW to reduce insanity of code
286	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
287
288	! included by RW for temporary comparison
289	real zdtnirco2(ngridmx,nlayermx) ! (K/s)
290
291	! included by RW to compute tracer column densities
292	real qcol(ngridmx,nqmx)
293
294	! included by RW for H2O precipitation
295	real zdtrain(ngridmx,nlayermx)
296	real zdqrain(ngridmx,nlayermx,nqmx)
297	real zdqsrain(ngridmx)
298	real zdqssnow(ngridmx)
299	real rainout(ngridmx)
300
301	! included by RW for H2O Manabe scheme
302	real dtmoist(ngridmx,nlayermx)
303	real dqmoist(ngridmx,nlayermx,nqmx)
304
305	real qvap(ngridmx,nlayermx)
306	real dqvaplscale(ngridmx,nlayermx)
307	real dqcldlscale(ngridmx,nlayermx)
308	real rneb_man(ngridmx,nlayermx)
309	real rneb_lsc(ngridmx,nlayermx)
310
311	! included by RW to account for surface cooling by evaporation
312	real dtsurfh2olat(ngridmx)
313
314	! included by RW to test energy conservation
315	real dEtot, dEtots, masse, vabs, dvabs
316	real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
317
318	real dItot, dVtot
319
320	! included by RW to test water conservation
321	real h2otot
322
323
324
325	! included by BC for evaporation
326	real qevap(ngridmx,nlayermx,nqmx)
327	real tevap(ngridmx,nlayermx)
328	real dqevap(ngridmx,nlayermx)
329	real dtevap(ngridmx,nlayermx)
330
331	! included by BC for hydrology
332	real hice(ngridmx)
333
334	! included by RW to test water conservation (by routine)
335	real dWtot, dWtots
336	real h2o_surf_all
337	logical watertest
338	save watertest
339
340	! included by RW for RH diagnostic
341	real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx)
342
343	! included by RW for hydrology
344	real dqs_hyd(ngridmx,nqmx)
345	real zdtsurf_hyd(ngridmx)
346
347	! included by RW for water cycle conservation tests
348	real icesrf,liqsrf,icecol,vapcol
349
350	! included by BC for double radiative transfer call
351	logical clearsky
352	real zdtsw1(ngridmx,nlayermx)
353	real zdtlw1(ngridmx,nlayermx)
354	real fluxsurf_lw1(ngridmx)
355	real fluxsurf_sw1(ngridmx)
356	real fluxtop_lw1(ngridmx)
357	real fluxabs_sw1(ngridmx)
358	real tau_col1(ngrid)
359	real OLR_nu1(ngrid,L_NSPECTI)
360	real OSR_nu1(ngrid,L_NSPECTV)
361	real tf, ntf
362
363	! included by BC for cloud fraction computations
364	real cloudfrac(ngridmx,nlayermx)
365	real totcloudfrac(ngridmx)
366
367	! included by RW for Tmax test
368	integer iTmax
369
370	! included by RW for vdifc water conservation test
371	real nconsMAX
372	real vdifcncons(ngridmx)
373	real cadjncons(ngridmx)
374	real qzero1D
375	save qzero1D
376
377	! double precision qsurf_hist(ngridmx,nqmx)
378	real qsurf_hist(ngridmx,nqmx)
379	save qsurf_hist
380
381	! included by RW for temp convadj conservation test
382	real playtest(nlayermx)
383	real plevtest(nlayermx)
384	real ttest(nlayermx)
385	real qtest(nlayermx)
386	integer igtest
387
388	! included by RW for runway greenhouse 1D study
389	real muvar(ngridmx,nlayermx+1)
390
391	! included by RW for variable H2O particle sizes
392	real reffH2O(ngridmx,nlayermx)
393	real reffcol(ngridmx,2)
394
395	! included by RW for flexible 3D energy conservation testing
396	real vdifcdE(ngridmx), madjdE(ngridmx), lscaledE(ngridmx)
397
398	! included by RW for sourceevol
399	real ice_initial(ngridmx)!, isoil
400	real delta_ice,ice_tot
401	real ice_min(ngridmx)
402	save ice_min
403	save ice_initial
404
405	integer num_run
406	logical ice_update
407	save ice_update
408
409	!=======================================================================
410
411
412	! 1. Initialisation
413	! -----------------
414
415	! 1.1 Initialisation only at first call
416	! ---------------------------------------
417	if (firstcall) then
418
419
420	! variables set to 0
421	! ~~~~~~~~~~~~~~~~~~
422	dtrad(:,:) = 0.0
423	fluxrad(:) = 0.0
424	tau_col(:) = 0.0
425	zdtsw(:,:) = 0.0
426	zdtlw(:,:) = 0.0
427
428	! initialize tracer names, indexes and properties
429	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
430	tracerdyn=tracer
431	if (tracer) then
432	call initracer()
433	endif ! end tracer
434
435	! read startfi (initial parameters)
436	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437	call phyetat0("startfi.nc",0,0,nsoilmx,nq, &
438	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
439	cloudfrac,totcloudfrac,hice)
440
441	if (pday.ne.day_ini) then
442	write(,) "ERROR in physiq.F90:"
443	write(,) "bad synchronization between physics and dynamics"
444	write(,) "dynamics day: ",pday
445	write(,) "physics day: ",day_ini
446	stop
447	endif
448
449	write (,) 'In physiq day_ini =', day_ini
450
451	! Initialize albedo and orbital calculation
452	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453	call surfini(ngrid,qsurf,albedo0)
454	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
455
456	do ig=1,ngrid
457	albedo(ig)=albedo0(ig)
458	enddo
459
460	if(tlocked)then
461	print*,'Planet is tidally locked at resonance n=',nres
462	print*,'Make sure you have the right rotation rate!!!'
463	endif
464
465	! initialize soil
466	! ~~~~~~~~~~~~~~~
467	if (callsoil) then
468	call soil(ngrid,nsoilmx,firstcall,inertiedat, &
469	ptimestep,tsurf,tsoil,capcal,fluxgrd)
470	else
471	print*,'WARNING! Thermal conduction in the soil turned off'
472	do ig=1,ngrid
473	capcal(ig)=1.e6
474	fluxgrd(ig)=0.
475	if(noradsurf)then
476	fluxgrd(ig)=10.0
477	endif
478	enddo
479	print*,'Flux from ground = ',fluxgrd,' W m^-2'
480	endif
481	icount=1
482
483	! decide whether to update ice at end of run
484	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485	ice_update=.false.
486	if(sourceevol)then
487	open(128,file='num_run',form='formatted')
488	read(128,*) num_run
489	close(128)
490
491	if(num_run.ne.0.and.mod(num_run,2).eq.0)then
492	!if(num_run.ne.0.and.mod(num_run,3).eq.0)then
493	print*,'Updating ice at end of this year!'
494	ice_update=.true.
495	ice_min(:)=1.e4
496	endif
497	endif
498
499	! define surface as continent or ocean
500	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
501	do ig=1,ngridmx
502	rnat(ig)=1
503
504	! if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
505	if(inertiedat(ig,1).gt.1.E4)then
506	rnat(ig)=0
507	endif
508	enddo
509
510	print*,'WARNING! Surface type currently decided by surface inertia'
511	print*,'This should be improved e.g. in newstart.F'
512
513
514	! initialise surface history variable
515	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
516	do ig=1,ngridmx
517	do iq=1,nqmx
518	qsurf_hist(ig,iq)=qsurf(ig,iq)
519	enddo
520	enddo
521
522	! initialise variable for dynamical heating diagnostic
523	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
524	ztprevious(:,:)=pt(:,:)
525
526	! Set temperature just above condensation temperature (for Early Mars)
527	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
528	if(nearco2cond) then
529	write(,)' WARNING! Starting at Tcond+1K'
530	do l=1, nlayer
531	do ig=1,ngrid
532	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
533	-pt(ig,l)) / ptimestep
534	enddo
535	enddo
536	endif
537
538	if(meanOLR)then
539	! to record global radiative balance
540	call system('rm -f rad_bal.out')
541	! to record global mean/max/min temperatures
542	call system('rm -f tem_bal.out')
543	! to record global hydrological balance
544	call system('rm -f h2o_bal.out')
545	endif
546
547	watertest=.false.
548	if(water)then
549	! initialise variables for water cycle
550
551	qzero1D=0.0
552
553	if(enertest)then
554	watertest = .true.
555	endif
556
557	if(ngrid.eq.1)then
558	qzero1D = 0.0
559	qsurf(1,igcm_h2o_vap) = qzero1D
560	do l=1, nlayer
561	pq(1,l,igcm_h2o_vap)=0.0
562	pq(1,l,igcm_h2o_ice)=0.0
563	enddo
564	endif
565
566	do ig=1,ngrid
567	qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
568	if(ice_update)then
569	ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
570	endif
571	enddo
572
573	endif
574	call su_watercycle ! even if we don't have a water cycle, we might
575	! need epsi for the wvp definitions in callcorrk.F
576
577	endif ! (end of "if firstcall")
578
579	! ---------------------------------------------------
580	! 1.2 Initializations done at every physical timestep:
581	! ---------------------------------------------------
582
583	if (ngrid.NE.ngridmx) then
584	print*,'STOP in PHYSIQ'
585	print*,'Probleme de dimensions :'
586	print*,'ngrid = ',ngrid
587	print*,'ngridmx = ',ngridmx
588	stop
589	endif
590
591	! Initialize various variables
592	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
593
594	pdu(:,:) = 0.0
595	pdv(:,:) = 0.0
596	! if ( (.not.nearco2cond).and.(.not.firstcall) ) then
597	if ( .not.nearco2cond ) then
598	pdt(:,:) = 0.0
599	endif ! this was the source of an evil bug...
600	pdq(:,:,:) = 0.0
601	pdpsrf(:) = 0.0
602	zflubid(:) = 0.0
603	zdtsurf(:) = 0.0
604	dqsurf(:,:) = 0.0
605
606	zday=pday+ptime ! compute time, in sols (and fraction thereof)
607
608	! Compute Stellar Longitude (Ls)
609	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
610	if (season) then
611	call stellarlong(zday,zls)
612	else
613	call stellarlong(float(day_ini),zls)
614	end if
615
616	! Compute geopotential between layers
617	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
618
619	do l=1,nlayer
620	do ig=1,ngrid
621	zzlay(ig,l)=pphi(ig,l)/g
622	enddo
623	enddo
624	do ig=1,ngrid
625	zzlev(ig,1)=0.
626	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
627	enddo
628	do l=2,nlayer
629	do ig=1,ngrid
630	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
631	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
632	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
633	enddo
634	enddo
635	! Potential temperature calculation may not be the same in physiq and dynamic...
636
637	! Compute potential temperature
638	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
639
640
641
642	do l=1,nlayer
643	do ig=1,ngrid
644	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
645	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
646	enddo
647	enddo
648
649	!-----------------------------------------------------------------------
650	! 2. Compute radiative tendencies
651	!-----------------------------------------------------------------------
652
653	if (callrad) then
654	if( mod(icount-1,iradia).eq.0.or.lastcall) then
655
656	! Compute local stellar zenith angles
657	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
658	call orbite(zls,dist_star,declin)
659
660	if (tlocked) then
661	ztim1=SIN(declin)
662	ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
663	ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
664
665	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
666	ztim1,ztim2,ztim3,mu0,fract)
667
668	elseif (diurnal) then
669	ztim1=SIN(declin)
670	ztim2=COS(declin)COS(2.pi*(zday-.5))
671	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
672
673	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
674	ztim1,ztim2,ztim3,mu0,fract)
675
676	else
677
678	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
679	! WARNING: this function appears not to work in 1D
680
681	endif
682
683	if (corrk) then
684
685	! a) Call correlated-k radiative transfer scheme
686	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
687
688	if(kastprof)then
689	print*,'kastprof should not = true here'
690	call abort
691	endif
692	muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
693
694	! standard callcorrk
695	clearsky=.false.
696	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
697	albedo,emis,mu0,pplev,pplay,pt, &
698	tsurf,fract,dist_star,aerosol,muvar, &
699	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
700	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
701	reffrad,tau_col,cloudfrac,totcloudfrac, &
702	clearsky,firstcall,lastcall)
703
704	! Option to call scheme once more for clear regions
705	if(CLFvarying)then
706
707	!!! ---> PROBLEMS WITH ALLOCATED ARRAYS
708	!!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
709	clearsky=.true.
710	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
711	albedo,emis,mu0,pplev,pplay,pt, &
712	tsurf,fract,dist_star,aerosol,muvar, &
713	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
714	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, &
715	reffrad,tau_col1,cloudfrac,totcloudfrac, &
716	clearsky,firstcall,lastcall)
717	clearsky = .false. !! just in case.
718
719	! Sum the fluxes and heating rates from cloudy/clear cases
720	do ig=1,ngrid
721	tf=totcloudfrac(ig)
722	ntf=1.-tf
723
724	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
725	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
726	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
727	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
728	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
729
730	do l=1,nlayer
731	zdtlw(ig,l) = ntfzdtlw1(ig,l) + tfzdtlw(ig,l)
732	zdtsw(ig,l) = ntfzdtsw1(ig,l) + tfzdtsw(ig,l)
733	enddo
734
735	do nw=1,L_NSPECTV
736	OSR_nu(ig,nw) = ntfOSR_nu1(ig,nw) + tfOSR_nu(ig,nw)
737	enddo
738	do nw=1,L_NSPECTI
739	OLR_nu(ig,nw) = ntfOLR_nu1(ig,nw) + tfOLR_nu(ig,nw)
740	enddo
741
742	enddo
743
744	endif !CLFvarying
745
746	! Radiative flux from the sky absorbed by the surface (W.m-2)
747	GSR=0.0
748	do ig=1,ngrid
749	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig) &
750	+fluxsurf_sw(ig)*(1.-albedo(ig))
751
752	if(noradsurf)then ! no lower surface; SW flux just disappears
753	GSR = GSR + fluxsurf_sw(ig)*area(ig)
754	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
755	endif
756
757	enddo
758	if(noradsurf)then
759	print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
760	endif
761
762	! Net atmospheric radiative heating rate (K.s-1)
763	do l=1,nlayer
764	do ig=1,ngrid
765	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
766	enddo
767	enddo
768
769
770	elseif(newtonian)then
771
772	! b) Call Newtonian cooling scheme
773	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
774	call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
775
776	do ig=1,ngrid
777	zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
778	enddo
779	! e.g. surface becomes proxy for 1st atmospheric layer ?
780
781	else
782
783	! c) Atmosphere has no radiative effect
784	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
785	do ig=1,ngrid
786	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_star**2
787	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
788	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
789	endif
790	fluxsurf_sw(ig) = fluxtop_dn(ig)
791	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
792	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
793	enddo ! radiation skips the atmosphere entirely
794
795	do l=1,nlayer
796	do ig=1,ngrid
797	dtrad(ig,l)=0.0
798	enddo
799	enddo ! hence no atmospheric radiative heating
800
801	endif ! if corrk
802
803	endif ! of if(mod(icount-1,iradia).eq.0)
804
805
806	! Transformation of the radiative tendencies
807	! ------------------------------------------
808
809	do ig=1,ngrid
810	zplanck(ig)=tsurf(ig)*tsurf(ig)
811	zplanck(ig)=emis(ig)stephanzplanck(ig)*zplanck(ig)
812	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
813	enddo
814
815	do l=1,nlayer
816	do ig=1,ngrid
817	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
818	enddo
819	enddo
820
821	!-------------------------
822	! test energy conservation
823	if(enertest)then
824	dEtotSW = 0.0
825	dEtotLW = 0.0
826	dEtotsSW = 0.0
827	dEtotsLW = 0.0
828	do ig = 1, ngrid
829	do l = 1, nlayer
830	masse = (pplev(ig,l) - pplev(ig,l+1))/g
831	dEtotSW = dEtotSW + cppmassezdtsw(ig,l)*area(ig)
832	dEtotLW = dEtotLW + cppmassezdtlw(ig,l)*area(ig)
833	enddo
834	dEtotsSW = dEtotsSW + fluxsurf_sw(ig)(1.-albedo(ig))area(ig)
835	dEtotsLW = dEtotsLW + emis(ig)fluxsurf_lw(ig)area(ig)-zplanck(ig)*area(ig)
836	enddo
837	dEtotSW = dEtotSW/totarea
838	dEtotLW = dEtotLW/totarea
839	dEtotsSW = dEtotsSW/totarea
840	dEtotsLW = dEtotsLW/totarea
841	print*,'---------------------------------------------------------------'
842	print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2'
843	print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2'
844	print*,'atmospheric energy change (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
845	print*,'In corrk SW surface energy change =',dEtotsSW,' W m-2'
846	print*,'In corrk LW surface energy change =',dEtotsLW,' W m-2'
847	print*,'surface energy change (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
848	endif
849	!-------------------------
850
851	endif ! of if (callrad)
852
853	!-----------------------------------------------------------------------
854	! 4. Vertical diffusion (turbulent mixing):
855	! -----------------------------------------
856
857	if (calldifv) then
858
859	do ig=1,ngrid
860	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
861	enddo
862
863	zdum1(:,:)=0.0
864	zdum2(:,:)=0.0
865	do l=1,nlayer
866	do ig=1,ngrid
867	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
868	enddo
869	enddo
870
871	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
872	ptimestep,capcal,lwrite, &
873	pplay,pplev,zzlay,zzlev,z0, &
874	pu,pv,zh,pq,tsurf,emis,qsurf, &
875	zdum1,zdum2,zdh,pdq,zflubid, &
876	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
877	zdqdif,zdqsdif,lastcall)
878
879	do l=1,nlayer
880	do ig=1,ngrid
881	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
882	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
883	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
884	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
885	enddo
886	enddo
887
888	do ig=1,ngrid
889	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
890	enddo
891
892	if (tracer) then
893	do iq=1, nq
894	do l=1,nlayer
895	do ig=1,ngrid
896	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
897	enddo
898	enddo
899	enddo
900	do iq=1, nq
901	do ig=1,ngrid
902	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
903	enddo
904	enddo
905
906	end if ! of if (tracer)
907
908	!-------------------------
909	! test energy conservation
910	if(enertest)then
911	dEtot=0.0
912	dEtots=0.0
913
914	vdifcdE(:)=0.0
915	do ig = 1, ngrid
916	do l = 1, nlayer
917	masse = (pplev(ig,l) - pplev(ig,l+1))/g
918	dEtot = dEtot + cppmassezdtdif(ig,l)*area(ig)
919
920	vabs = sqrt(pdu(ig,l)2 + pdv(ig,l)2)
921	dvabs = sqrt(zdudif(ig,l)2 + zdvdif(ig,l)2)
922	dEtot = dEtot + massevabsdvabs*area(ig)
923
924	vdifcdE(ig) = vdifcdE(ig) + massevabsdvabs
925
926	enddo
927	dEtot = dEtot - zflubid(ig)*area(ig) ! subtract flux from ground
928
929	dEtots = dEtots + capcal(ig)zdtsdif(ig)area(ig)
930	vdifcdE(ig) = vdifcdE(ig) + capcal(ig)*zdtsdif(ig)
931	enddo
932	dEtot = dEtot/totarea
933	dEtots = dEtots/totarea
934	print*,'In difv atmospheric energy change =',dEtot,' W m-2'
935	print*,'In difv surface energy change =',dEtots,' W m-2'
936	!print*,'Note we ignore the wind change...'
937	! and latent heat for that matter...
938	endif
939	!-------------------------
940
941
942	!-------------------------
943	! test water conservation
944	if(watertest.and.water)then
945	dWtot=0.0
946	dWtots=0.0
947	nconsMAX=0.0
948	do ig = 1, ngrid
949
950	vdifcncons(ig)=0.0
951	do l = 1, nlayer
952	masse = (pplev(ig,l) - pplev(ig,l+1))/g
953
954	iq = igcm_h2o_vap
955	dWtot = dWtot + massezdqdif(ig,l,iq)ptimestep*area(ig)
956	vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
957
958	iq = igcm_h2o_ice
959	dWtot = dWtot + massezdqdif(ig,l,iq)ptimestep*area(ig)
960	vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
961
962	enddo
963
964	iq = igcm_h2o_vap
965	dWtots = dWtots + zdqsdif(ig,iq)ptimesteparea(ig)
966	vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
967
968	iq = igcm_h2o_ice
969	dWtots = dWtots + zdqsdif(ig,iq)ptimesteparea(ig)
970	vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
971
972	if(vdifcncons(ig).gt.nconsMAX)then
973	nconsMAX=vdifcncons(ig)
974	endif
975
976	enddo
977
978	dWtot = dWtot/totarea
979	dWtots = dWtots/totarea
980	print*,'---------------------------------------------------------------'
981	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
982	print*,'In difv surface water change =',dWtots,' kg m-2'
983	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
984	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
985
986
987	endif
988	!-------------------------
989
990	else
991
992	if(.not.newtonian)then
993
994	do ig=1,ngrid
995	zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
996	enddo
997
998	endif
999
1000	endif ! of if (calldifv)
1001
1002
1003	!-----------------------------------------------------------------------
1004	! 5. Dry convective adjustment:
1005	! -----------------------------
1006
1007	if(calladj) then
1008
1009	do l=1,nlayer
1010	do ig=1,ngrid
1011	zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
1012	enddo
1013	enddo
1014	zduadj(:,:)=0.0
1015	zdvadj(:,:)=0.0
1016	zdhadj(:,:)=0.0
1017
1018
1019	call convadj(ngrid,nlayer,nq,ptimestep, &
1020	pplay,pplev,zpopsk, &
1021	pu,pv,zh,pq, &
1022	pdu,pdv,zdh,pdq, &
1023	zduadj,zdvadj,zdhadj, &
1024	zdqadj)
1025
1026	do l=1,nlayer
1027	do ig=1,ngrid
1028	pdu(ig,l) = pdu(ig,l) + zduadj(ig,l)
1029	pdv(ig,l) = pdv(ig,l) + zdvadj(ig,l)
1030	pdt(ig,l) = pdt(ig,l) + zdhadj(ig,l)*zpopsk(ig,l)
1031	zdtadj(ig,l) = zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1032	enddo
1033	enddo
1034
1035	if(tracer) then
1036	do iq=1, nq
1037	do l=1,nlayer
1038	do ig=1,ngrid
1039	pdq(ig,l,iq) = pdq(ig,l,iq) + zdqadj(ig,l,iq)
1040	enddo
1041	enddo
1042	enddo
1043	end if
1044
1045	!-------------------------
1046	! test energy conservation
1047	if(enertest)then
1048	dEtot=0.0
1049	do ig = 1, ngrid
1050	do l = 1, nlayer
1051	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1052	dEtot = dEtot + cppmassezdtadj(ig,l)*area(ig)
1053	enddo
1054	enddo
1055	dEtot=dEtot/totarea
1056	print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1057	endif
1058	!-------------------------
1059
1060	!-------------------------
1061	! test water conservation
1062	if(watertest)then
1063	dWtot=0.0
1064	do ig = 1, ngrid
1065	cadjncons(ig)=0.0
1066	do l = 1, nlayer
1067	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1068
1069	iq = igcm_h2o_vap
1070	dWtot = dWtot + massezdqadj(ig,l,iq)ptimestep*area(ig)
1071	cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
1072
1073	iq = igcm_h2o_ice
1074	dWtot = dWtot + massezdqadj(ig,l,iq)ptimestep*area(ig)
1075	cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
1076
1077	enddo
1078
1079	enddo
1080	dWtot=dWtot/totarea
1081	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
1082	endif
1083	!-------------------------
1084
1085	endif ! of if(calladj)
1086
1087	!-----------------------------------------------------------------------
1088	! 6. Carbon dioxide condensation-sublimation:
1089	! -------------------------------------------
1090
1091	if (co2cond) then
1092	if (.not.tracer) then
1093	print*,'We need a CO2 ice tracer to condense CO2'
1094	call abort
1095	endif
1096
1097	call condense_cloud(ngrid,nlayer,nq,ptimestep, &
1098	capcal,pplay,pplev,tsurf,pt, &
1099	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1100	qsurf(1,igcm_co2_ice),albedo,emis, &
1101	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1102	zdqc,reffrad)
1103
1104	do l=1,nlayer
1105	do ig=1,ngrid
1106	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1107	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1108	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1109	enddo
1110	enddo
1111	do ig=1,ngrid
1112	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1113	enddo
1114
1115	do iq=1,nq ! should use new notation here !
1116	do l=1,nlayer
1117	do ig=1,ngrid
1118	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1119	enddo
1120	enddo
1121	enddo
1122	! Note: we do not add surface co2ice tendency
1123	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
1124
1125	!-------------------------
1126	! test energy conservation
1127	if(enertest)then
1128	dEtot=0.0
1129	dEtots=0.0
1130	do ig = 1, ngrid
1131	do l = 1, nlayer
1132	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1133	dEtot = dEtot + cppmassezdtc(ig,l)*area(ig)
1134	enddo
1135	dEtots = dEtots + capcal(ig)zdtsurfc(ig)area(ig)
1136	enddo
1137	dEtot=dEtot/totarea
1138	dEtots=dEtots/totarea
1139	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
1140	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
1141	endif
1142	!-------------------------
1143
1144	endif ! of if (co2cond)
1145
1146
1147	!-----------------------------------------------------------------------
1148	! 7. Specific parameterizations for tracers
1149	! -----------------------------------------
1150
1151	if (tracer) then
1152
1153	! 7a. Water and ice
1154	! ---------------
1155	if (water) then
1156
1157	! ----------------------------------------
1158	! Water ice condensation in the atmosphere
1159	! ----------------------------------------
1160	if(watercond)then
1161
1162	if(RLVTT.gt.1.e-8)then
1163
1164	! Re-evaporate cloud water/ice
1165	call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1166	DO l = 1, nlayer
1167	DO ig = 1, ngrid
1168	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1169	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1170	pdt(ig,l) = pdt(ig,l)+dtevap(ig,l)
1171	enddo
1172	enddo
1173
1174	call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1175	do l=1,nlayer
1176	do ig=1,ngrid
1177	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
1178	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
1179	pdt(ig,l) = pdt(ig,l)+dtmoist(ig,l)
1180	enddo
1181	enddo
1182
1183	!-------------------------
1184	! test energy conservation
1185	if(enertest)then
1186	dEtot=0.0
1187	madjdE(:)=0.0
1188	do ig = 1, ngrid
1189	do l = 1, nlayer
1190	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1191	dEtot = dEtot + cppmasse(dtmoist(ig,l)+dtevap(ig,l))*area(ig)
1192	madjdE(ig) = madjdE(ig) + cppmasse(dtmoist(ig,l)+dtevap(ig,l))
1193	enddo
1194	enddo
1195	dEtot=dEtot/totarea
1196	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
1197	endif
1198	!-------------------------
1199
1200	!-------------------------
1201	! test water conservation
1202	if(watertest)then
1203	dWtot=0.0
1204	do ig = 1, ngrid
1205	do iq = 1 , nq
1206	do l = 1, nlayer
1207	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1208	dWtot = dWtot + massedqmoist(ig,l,igcm_h2o_vap)area(ig)*ptimestep
1209	dWtot = dWtot + massedqmoist(ig,l,igcm_h2o_ice)area(ig)*ptimestep
1210	enddo
1211	enddo
1212	enddo
1213	dWtot=dWtot/totarea
1214	print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
1215	endif
1216	!-------------------------
1217
1218
1219	endif
1220
1221
1222	! Re-evaporate cloud water/ice
1223	call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1224	do l = 1, nlayer
1225	do ig = 1, ngrid
1226	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1227	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1228	pdt(ig,l) = pdt(ig,l)+dtevap(ig,l)
1229	enddo
1230	enddo ! note: we use qevap but not tevap in largescale/moistadj
1231	! otherwise is a big mess
1232
1233	call largescale(ptimestep,pplev,pplay,pt,qevap, & ! a bug was here!
1234	pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
1235	do l=1,nlayer
1236	do ig=1,ngrid
1237	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
1238	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
1239	pdt(ig,l) = pdt(ig,l)+dtlscale(ig,l)
1240
1241	if(.not.aerofixed)then
1242	reffrad(ig,l,2)=reffH2O(ig,l)
1243	endif
1244
1245	enddo
1246	enddo
1247
1248	!-------------------------
1249	! test energy conservation
1250	if(enertest)then
1251	dEtot=0.0
1252	lscaledE(:)=0.0
1253	do ig = 1, ngrid
1254	do l = 1, nlayer
1255	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1256	dEtot = dEtot + cppmasse(dtlscale(ig,l)+dtevap(ig,l))*area(ig)
1257	lscaledE(ig) = lscaledE(ig) + cppmasse(dtlscale(ig,l)+dtevap(ig,l))
1258	enddo
1259	enddo
1260	dEtot=dEtot/totarea
1261	print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
1262	endif
1263	!-------------------------
1264
1265	!-------------------------
1266	! test water conservation
1267	if(watertest)then
1268	dWtot=0.0
1269	do ig = 1, ngrid
1270	do iq = 1 , nq
1271	do l = 1, nlayer
1272	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1273	dWtot = dWtot + massedqvaplscale(ig,l)area(ig)*ptimestep
1274	dWtot = dWtot + massedqcldlscale(ig,l)area(ig)*ptimestep
1275	enddo
1276	enddo
1277	enddo
1278	dWtot=dWtot/totarea
1279	print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
1280	endif
1281	!-------------------------
1282
1283	! compute cloud fraction
1284	do l = 1, nlayer
1285	do ig = 1,ngrid
1286	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1287	enddo
1288	enddo
1289
1290	! compute total cloud fraction in column
1291	call totalcloudfrac(cloudfrac,totcloudfrac)
1292
1293	endif ! of if (watercondense)
1294
1295
1296	! --------------------------------
1297	! Water ice / liquid precipitation
1298	! --------------------------------
1299	if(waterrain)then
1300
1301	zdqrain(:,:,:) = 0.0
1302	zdqsrain(:) = 0.0
1303	zdqssnow(:) = 0.0
1304
1305	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
1306	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
1307
1308	do l=1,nlayer
1309	do ig=1,ngrid
1310	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
1311	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
1312	pdt(ig,l) = pdt(ig,l)+zdtrain(ig,l)
1313	enddo
1314	enddo
1315
1316	do ig=1,ngrid
1317	dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
1318	dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
1319	rainout(ig) = zdqsrain(ig)+zdqssnow(ig) ! diagnostic
1320	enddo
1321
1322	!-------------------------
1323	! test energy conservation
1324	if(enertest)then
1325	dEtot=0.0
1326	do ig = 1, ngrid
1327	do l = 1, nlayer
1328	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1329	dEtot = dEtot + cppmassezdtrain(ig,l)*area(ig)
1330	enddo
1331	enddo
1332	dEtot=dEtot/totarea
1333	print*,'In rain atmospheric energy change =',dEtot,' W m-2'
1334	endif
1335	!-------------------------
1336
1337
1338	!-------------------------
1339	! test energy conservation
1340	if(enertest)then
1341	dEtot=0.0
1342	do ig = 1, ngrid
1343	do l = 1, nlayer
1344	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1345	dEtot = dEtot + cppmassezdtrain(ig,l)*area(ig)
1346	enddo
1347	enddo
1348	dEtot=dEtot/totarea
1349	print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
1350
1351	dEtot=0.0
1352	do ig = 1, ngrid
1353	do l = 1, nlayer
1354	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1355	dItot = dItot + massezdqrain(ig,l,igcm_h2o_ice)area(ig)
1356	dVtot = dVtot + massezdqrain(ig,l,igcm_h2o_vap)area(ig)
1357	enddo
1358	dItot = dItot + zdqssnow(ig)*area(ig)
1359	dVtot = dVtot + zdqsrain(ig)*area(ig)
1360	enddo
1361	dEtot=(dItotRLVTT/cpp + dVtotRLVTT/cpp)/totarea
1362	print,'In rain dItot =',dItotRLVTT/(cpp*totarea),' W m-2'
1363	print,'In rain dVtot =',dVtotRLVTT/(cpp*totarea),' W m-2'
1364	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
1365	endif
1366	!-------------------------
1367
1368	!-------------------------
1369	! test water conservation
1370	if(watertest)then
1371	dWtot=0.0
1372	dWtots=0.0
1373	do ig = 1, ngrid
1374	!do iq = 1 , nq
1375	do l = 1, nlayer
1376	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1377	dWtot = dWtot + massezdqrain(ig,l,igcm_h2o_vap)area(ig)*ptimestep
1378	dWtot = dWtot + massezdqrain(ig,l,igcm_h2o_ice)area(ig)*ptimestep
1379	enddo
1380	!enddo
1381	dWtots = dWtots + (zdqsrain(ig)+zdqssnow(ig))area(ig)ptimestep
1382	enddo
1383	dWtot=dWtot/totarea
1384	dWtots=dWtots/totarea
1385	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
1386	print*,'In rain surface water change =',dWtots,' kg m-2'
1387	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
1388	endif
1389	!-------------------------
1390
1391	end if ! of if (waterrain)
1392	end if ! of if (water)
1393
1394
1395	! 7c. Aerosol particles
1396	! -------------------
1397	! -------------
1398	! Sedimentation
1399	! -------------
1400	if (sedimentation) then
1401	zdqsed(:,:,:) = 0.0
1402	zdqssed(:,:) = 0.0
1403
1404
1405	!-------------------------
1406	! find qtot
1407	if(watertest)then
1408	dWtot=0.0
1409	dWtots=0.0
1410	iq=3
1411	do ig = 1, ngrid
1412	do l = 1, nlayer
1413	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1414	dWtot = dWtot + massepq(ig,l,iq)area(ig)*ptimestep
1415	dWtots = dWtots + massepdq(ig,l,iq)area(ig)*ptimestep
1416	enddo
1417	enddo
1418	dWtot=dWtot/totarea
1419	dWtots=dWtots/totarea
1420	print*,'Before sedim pq =',dWtot,' kg m-2'
1421	print*,'Before sedim pdq =',dWtots,' kg m-2'
1422	endif
1423	!-------------------------
1424
1425	call callsedim(ngrid,nlayer,ptimestep, &
1426	pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
1427
1428	!-------------------------
1429	! find qtot
1430	if(watertest)then
1431	dWtot=0.0
1432	dWtots=0.0
1433	iq=3
1434	do ig = 1, ngrid
1435	do l = 1, nlayer
1436	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1437	dWtot = dWtot + massepq(ig,l,iq)area(ig)*ptimestep
1438	dWtots = dWtots + massepdq(ig,l,iq)area(ig)*ptimestep
1439	enddo
1440	enddo
1441	dWtot=dWtot/totarea
1442	dWtots=dWtots/totarea
1443	print*,'After sedim pq =',dWtot,' kg m-2'
1444	print*,'After sedim pdq =',dWtots,' kg m-2'
1445	endif
1446	!-------------------------
1447
1448	do iq=1,nq
1449	! for now, we only allow H2O ice to sediment
1450	! and as in rain.F90, whether it falls as rain or snow depends
1451	! only on the surface temperature
1452	do ig=1,ngrid
1453	do l=1,nlayer
1454	pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
1455	enddo
1456	dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
1457	enddo
1458	enddo
1459
1460	!-------------------------
1461	! test water conservation
1462	if(watertest)then
1463	dWtot=0.0
1464	dWtots=0.0
1465	do iq=1,nq
1466	do ig = 1, ngrid
1467	do l = 1, nlayer
1468	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1469	dWtot = dWtot + massezdqsed(ig,l,iq)area(ig)*ptimestep
1470	enddo
1471	dWtots = dWtots + zdqssed(ig,iq)area(ig)ptimestep
1472	enddo
1473	enddo
1474	dWtot=dWtot/totarea
1475	dWtots=dWtots/totarea
1476	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
1477	print*,'In sedim surface ice change =',dWtots,' kg m-2'
1478	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
1479	endif
1480	!-------------------------
1481
1482	end if ! of if (sedimentation)
1483
1484
1485	! 7d. Updates
1486	! ---------
1487
1488	! ---------------------------------
1489	! Updating tracer budget on surface
1490	! ---------------------------------
1491
1492	if(hydrology)then
1493
1494	call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
1495	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
1496	! note: for now, also changes albedo in the subroutine
1497
1498	do ig=1,ngrid
1499	zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
1500	do iq=1,nq
1501	qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
1502	enddo
1503	enddo
1504	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
1505
1506	!-------------------------
1507	! test energy conservation
1508	if(enertest)then
1509	dEtot=0.0
1510	do ig = 1, ngrid
1511	dEtots = dEtots + capcal(ig)zdtsurf_hyd(ig)area(ig)
1512	enddo
1513	dEtot=dEtot/totarea
1514	print*,'In hydrol atmospheric energy change =',dEtot,' W m-2'
1515	endif
1516	!-------------------------
1517
1518	!-------------------------
1519	! test water conservation
1520	if(watertest)then
1521	dWtots=0.0
1522	do ig = 1, ngrid
1523	dWtots = dWtots + dqs_hyd(ig,igcm_h2o_ice)area(ig)ptimestep
1524	enddo
1525	dWtots=dWtots/totarea
1526	print*,'In hydrol surface ice change =',dWtots,' kg m-2'
1527	dWtots=0.0
1528	do ig = 1, ngrid
1529	dWtots = dWtots + dqs_hyd(ig,igcm_h2o_vap)area(ig)ptimestep
1530	enddo
1531	dWtots=dWtots/totarea
1532	print*,'In hydrol surface water change =',dWtots,' kg m-2'
1533	print*,'---------------------------------------------------------------'
1534	endif
1535	!-------------------------
1536
1537	ELSE ! of if (hydrology)
1538
1539	do iq=1,nq
1540	do ig=1,ngrid
1541	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1542	enddo
1543	enddo
1544
1545	END IF ! of if (hydrology)
1546
1547	! Add qsurf to qsurf_hist, which is what we save in
1548	! diagfi.nc etc. At the same time, we set the water
1549	! content of ocean gridpoints back to zero, in order
1550	! to avoid rounding errors in vdifc, rain
1551	do ig = 1, ngrid
1552	do iq = 1, nq
1553	if(iq.eq.igcm_h2o_vap .and. rnat(ig).eq.0)then ! if liquid water and terrain = ocean
1554	qsurf_hist(ig,iq) = qsurf(ig,iq)
1555	!qsurf(ig,iq) = qcol(ig,iq)
1556	! the value of qsurf we choose here makes NO DIFFERENCE TO ANYTHING AT ALL
1557	else
1558	qsurf_hist(ig,iq) = qsurf(ig,iq)
1559	endif
1560	enddo
1561	enddo
1562
1563	if(ice_update)then
1564	do ig = 1, ngrid
1565	ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
1566	enddo
1567	endif
1568
1569	endif ! of if (tracer)
1570
1571	!-----------------------------------------------------------------------
1572	! 9. Surface and sub-surface soil temperature
1573	!-----------------------------------------------------------------------
1574
1575
1576	! Increment surface temperature
1577	do ig=1,ngrid
1578	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1579	enddo
1580
1581	! Compute soil temperatures and subsurface heat flux
1582	if (callsoil) then
1583	call soil(ngrid,nsoilmx,.false.,inertiedat, &
1584	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1585	endif
1586
1587	!-------------------------
1588	! test energy conservation
1589	if(enertest)then
1590	dEtots=0.0
1591	do ig = 1, ngrid
1592	dEtots = dEtots + capcal(ig)zdtsurf(ig)area(ig)
1593	enddo
1594	dEtots=dEtots/totarea
1595	print*,'Surface energy change=',dEtots,' W m-2'
1596	endif
1597	!-------------------------
1598
1599	!-----------------------------------------------------------------------
1600	! 10. Perform diagnostics and write output files
1601	!-----------------------------------------------------------------------
1602
1603	! -------------------------------
1604	! Dynamical fields incrementation
1605	! -------------------------------
1606	! For output only: the actual model integration is performed in the dynamics
1607
1608	! temperature, zonal and meridional wind
1609	do l=1,nlayer
1610	do ig=1,ngrid
1611	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
1612	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
1613	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
1614
1615	! diagnostic
1616	zdtdyn(ig,l) = ztprevious(ig,l)-pt(ig,l)
1617	ztprevious(ig,l) = zt(ig,l)
1618	enddo
1619	enddo
1620
1621	if(firstcall)then
1622	zdtdyn(:,:)=0.0
1623	endif
1624
1625	! dynamical heating diagnostic
1626	fluxdyn(:)=0.
1627	do ig=1,ngrid
1628	do l=1,nlayer
1629	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep) &
1630	(pplev(ig,l)-pplev(ig,l+1))cpp/g
1631	enddo
1632	enddo
1633
1634	! tracers
1635	do iq=1, nq
1636	do l=1,nlayer
1637	do ig=1,ngrid
1638	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
1639	enddo
1640	enddo
1641	enddo
1642
1643	! surface pressure
1644	do ig=1,ngrid
1645	ps(ig) = pplev(ig,1) + pdpsrf(ig)*ptimestep
1646	enddo
1647
1648	! pressure
1649	do l=1,nlayer
1650	do ig=1,ngrid
1651	zplev(ig,l) = pplev(ig,l)/pplev(ig,1)*ps(ig)
1652	zplay(ig,l) = pplay(ig,l)/pplev(ig,1)*ps(ig)
1653	enddo
1654	enddo
1655
1656	! ---------------------------------------------------------
1657	! Surface and soil temperature information
1658	! ---------------------------------------------------------
1659
1660	Ts1 = 0.0
1661	Ts2 = 99999.9
1662	Ts3 = 0.0
1663	TsS = 0.0 ! mean temperature at bottom soil layer
1664	do ig=1,ngrid
1665	Ts1 = Ts1 + area(ig)*tsurf(ig)
1666	Ts2 = min(Ts2,tsurf(ig))
1667	Ts3 = max(Ts3,tsurf(ig))
1668	TsS = TsS + area(ig)*tsoil(ig,nsoilmx)
1669	end do
1670	Ts1=Ts1/totarea
1671	TsS=TsS/totarea
1672	if(callsoil)then
1673	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1674	print*,Ts1,Ts2,Ts3,TsS
1675	else
1676	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1677	print*,Ts1,Ts2,Ts3
1678	endif
1679
1680	! ---------------------------------------------------------
1681	! Check the energy balance of the simulation during the run
1682	! ---------------------------------------------------------
1683
1684	if(corrk)then
1685
1686	ISR = 0.0
1687	ASR = 0.0
1688	OLR = 0.0
1689	GND = 0.0
1690	DYN = 0.0
1691	do ig=1,ngrid
1692	ISR = ISR + area(ig)*fluxtop_dn(ig)
1693	ASR = ASR + area(ig)*fluxabs_sw(ig)
1694	OLR = OLR + area(ig)*fluxtop_lw(ig)
1695	GND = GND + area(ig)*fluxgrd(ig)
1696	if(.not.callsoil) GND=GND+ area(ig)*fluxrad(ig)
1697	DYN = DYN + area(ig)*fluxdyn(ig)
1698
1699	if(fluxtop_dn(ig).lt.0.0)then
1700	print*,'fluxtop_dn has gone crazy'
1701	print*,'fluxtop_dn=',fluxtop_dn(ig)
1702	print*,'tau_col=',tau_col(ig)
1703	print*,'aerosol=',aerosol(ig,:,:)
1704	print*,'temp= ',pt(ig,:)
1705	print*,'pplay= ',pplay(ig,:)
1706	call abort
1707	endif
1708	end do
1709
1710	if(ngridmx.eq.1)then
1711	DYN=0.0
1712	endif
1713
1714	print*,' ISR ASR OLR GND DYN [W m^-2]'
1715	print*, ISR/totarea,ASR/totarea,OLR/totarea,GND/totarea,DYN/totarea
1716
1717	if(enertest)then
1718	print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
1719	print,'LW energy balance LW++ + ASR* = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
1720	print,'LW energy balance LW++ - OLR* = ',dEtotLW+dEtotsLW+OLR/totarea,' W m-2'
1721	endif
1722
1723	if(meanOLR)then
1724	if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
1725	! to record global radiative balance
1726	open(92,file="rad_bal.out",form='formatted',access='append')
1727	write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
1728	close(92)
1729	open(93,file="tem_bal.out",form='formatted',access='append')
1730	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1731	close(93)
1732	endif
1733	endif
1734
1735	endif
1736
1737	! ------------------------------------------------------------------
1738	! Diagnostic to test radiative-convective timescales in code
1739	! ------------------------------------------------------------------
1740	if(testradtimes)then
1741	open(38,file="tau_phys.out",form='formatted',access='append')
1742	ig=1
1743	do l=1,nlayer
1744	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1745	enddo
1746	close(38)
1747	print*,'As testradtimes enabled, exiting physics on first call'
1748	call abort
1749	endif
1750
1751	! ---------------------------------------------------------
1752	! Compute column amounts (kg m-2) if tracers are enabled
1753	! ---------------------------------------------------------
1754	if(tracer)then
1755	qcol(:,:)=0.0
1756	do iq=1,nq
1757	do ig=1,ngrid
1758	do l=1,nlayer
1759	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) * &
1760	(pplev(ig,l) - pplev(ig,l+1)) / g
1761	enddo
1762	enddo
1763	enddo
1764
1765	! not generalised for arbitrary aerosols yet!!!
1766	reffcol(:,:)=0.0
1767	do ig=1,ngrid
1768	do l=1,nlayer
1769	if(co2cond)then
1770	reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &
1771	reffrad(ig,l,1) * &
1772	(pplev(ig,l) - pplev(ig,l+1)) / g
1773	endif
1774	if(water)then
1775	reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &
1776	reffrad(ig,l,2) * &
1777	(pplev(ig,l) - pplev(ig,l+1)) / g
1778	endif
1779	enddo
1780	enddo
1781
1782	endif
1783
1784	! ---------------------------------------------------------
1785	! Test for water conservation if water is enabled
1786	! ---------------------------------------------------------
1787
1788	if(water)then
1789
1790	icesrf = 0.0
1791	liqsrf = 0.0
1792	icecol = 0.0
1793	vapcol = 0.0
1794
1795	h2otot = 0.0
1796	do ig=1,ngrid
1797
1798	icesrf = icesrf + area(ig)*qsurf_hist(ig,igcm_h2o_ice)
1799	liqsrf = liqsrf + area(ig)*qsurf_hist(ig,igcm_h2o_vap)
1800	icecol = icecol + area(ig)*qcol(ig,igcm_h2o_ice)
1801	vapcol = vapcol + area(ig)*qcol(ig,igcm_h2o_vap)
1802
1803	h2otot = h2otot + area(ig)* &
1804	(qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap) &
1805	+qsurf_hist(ig,igcm_h2o_ice)+qsurf_hist(ig,igcm_h2o_vap))
1806	end do
1807
1808	print*,' Total water amount [kg m^-2]: ',h2otot/totarea
1809	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
1810	print*, icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
1811
1812	if(meanOLR)then
1813	if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
1814	! to record global water balance
1815	open(98,file="h2o_bal.out",form='formatted',access='append')
1816	write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
1817	close(98)
1818	endif
1819	endif
1820
1821	endif
1822
1823	! ---------------------------------------------------------
1824	! Calculate RH for diagnostic if water is enabled
1825	! ---------------------------------------------------------
1826
1827	if(water)then
1828	do l = 1, nlayer
1829	do ig = 1, ngrid
1830	call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
1831	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
1832	enddo
1833	enddo
1834
1835	! compute maximum possible H2O column amount (100% saturation)
1836	do ig=1,ngrid
1837	H2Omaxcol(ig)=0.0
1838	do l=1,nlayer
1839	H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) * &
1840	(pplev(ig,l) - pplev(ig,l+1))/g
1841	enddo
1842	enddo
1843
1844	endif
1845
1846
1847	print*,''
1848	print,'--> Ls =',zls180./pi
1849	! -------------------------------------------------------------------
1850	! Writing NetCDF file "RESTARTFI" at the end of the run
1851	! -------------------------------------------------------------------
1852	! Note: 'restartfi' is stored just before dynamics are stored
1853	! in 'restart'. Between now and the writting of 'restart',
1854	! there will have been the itau=itau+1 instruction and
1855	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1856	! thus we store for time=time+dtvr
1857
1858	if(lastcall) then
1859	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1860
1861
1862	! Update surface ice distribution to iterate to steady state if requested
1863	if(ice_update)then
1864
1865	do ig = 1, ngrid
1866
1867	delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
1868
1869	! add multiple years of evolution
1870	qsurf_hist(ig,igcm_h2o_ice) = &
1871	!qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0
1872	qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
1873
1874	! if ice has gone -ve, set to zero
1875	if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
1876	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1877	!qsurf_hist(ig,igcm_h2o_vap) = 0.0
1878	endif
1879
1880	! if ice is seasonal, set to zero (NEW)
1881	if(ice_min(ig).lt.0.01)then
1882	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1883	!qsurf_hist(ig,igcm_h2o_vap) = 0.0
1884	endif
1885
1886	enddo
1887
1888	! enforce ice conservation
1889	ice_tot=0.0
1890	do ig = 1, ngrid
1891	ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
1892	enddo
1893	do ig = 1, ngrid
1894	qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
1895	enddo
1896
1897	endif
1898
1899	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1900	call physdem1("restartfi.nc",long,lati,nsoilmx,nq, &
1901	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1902	area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1903	cloudfrac,totcloudfrac,hice)
1904	endif
1905
1906	! -----------------------------------------------------------------
1907	! Saving statistics :
1908	! -----------------------------------------------------------------
1909	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1910	! which can later be used to make the statistic files of the run:
1911	! "stats") only possible in 3D runs !
1912
1913
1914	if (callstats) then
1915
1916	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1917	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1918	call wstats(ngrid,"fluxsurf_lw", &
1919	"Thermal IR radiative flux to surface","W.m-2",2, &
1920	fluxsurf_lw)
1921	! call wstats(ngrid,"fluxsurf_sw", &
1922	! "Solar radiative flux to surface","W.m-2",2, &
1923	! fluxsurf_sw_tot)
1924	call wstats(ngrid,"fluxtop_lw", &
1925	"Thermal IR radiative flux to space","W.m-2",2, &
1926	fluxtop_lw)
1927	! call wstats(ngrid,"fluxtop_sw", &
1928	! "Solar radiative flux to space","W.m-2",2, &
1929	! fluxtop_sw_tot)
1930
1931	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1932	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1933	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1934
1935	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1936	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1937	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1938	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1939	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1940
1941	if (tracer) then
1942	do iq=1,nq
1943	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1944	call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1945	'kg m^-2',2,qsurf(1,iq) )
1946
1947	call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
1948	'kg m^-2',2,qcol(1,iq) )
1949	call wstats(ngridmx,trim(noms(iq))//'_reff', &
1950	trim(noms(iq))//'_reff', &
1951	'm',3,reffrad(1,1,iq))
1952	end do
1953	if (water) then
1954	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
1955	call wstats(ngrid,"vmr_h2ovapor", &
1956	"H2O vapour volume mixing ratio","mol/mol", &
1957	3,vmr)
1958	endif ! of if (water)
1959
1960	endif !tracer
1961
1962	if(lastcall) then
1963	write (,) "Writing stats..."
1964	call mkstats(ierr)
1965	endif
1966	endif !if callstats
1967
1968
1969	! ----------------------------------------------------------------------
1970	! output in netcdf file "DIAGFI", containing any variable for diagnostic
1971	! (output with period "ecritphy", set in "run.def")
1972	! ----------------------------------------------------------------------
1973	! writediagfi can also be called from any other subroutine for any variable.
1974	! but its preferable to keep all the calls in one place...
1975
1976	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1977	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1978	call writediagfi(ngrid,"temp","temperature","K",3,zt)
1979	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1980	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1981	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1982	call writediagfi(ngrid,'p','Pressure','Pa',3,pplay)
1983
1984	! Total energy balance diagnostics
1985	if(callrad.and.(.not.newtonian))then
1986	call writediagfi(ngrid,'ALB','Surface albedo',' ',2,albedo)
1987	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1988	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1989	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1990	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
1991	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1992	!call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
1993	!call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
1994	!call writediagfi(ngrid,"vdifcdE","heat from vdifc surface","W m-2",2,vdifcdE)
1995	endif
1996
1997	! Temporary inclusions for heating diagnostics
1998	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1999	! call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2000	! call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2001	! call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2002
2003	! debugging
2004	!call writediagfi(ngrid,"vdifNC","H2O loss vdifc","kg m-2 s-1",2,vdifcncons)
2005	!call writediagfi(ngrid,"cadjNC","H2O loss convadj","kg m-2 s-1",2,cadjncons)
2006	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2007	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2008
2009	! Output aerosols
2010	if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
2011	if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
2012	if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
2013	if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
2014
2015	! Output tracers
2016	if (tracer) then
2017	do iq=1,nq
2018	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2019	! call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2020	! 'kg m^-2',2,qsurf(1,iq) )
2021	call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2022	'kg m^-2',2,qsurf_hist(1,iq) )
2023	call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2024	'kg m^-2',2,qcol(1,iq) )
2025
2026	if(water)then
2027	call writediagfi(ngridmx,"H2Omaxcol","max. poss. H2O column","kg m^-2",2,H2Omaxcol)
2028	endif
2029
2030	if(watercond)then
2031	!call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2032	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2033	endif
2034
2035	if(waterrain)then
2036	call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
2037	call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
2038	endif
2039
2040	if(hydrology)then
2041	call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
2042	endif
2043
2044	if(ice_update)then
2045	call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
2046	call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
2047	endif
2048
2049	do ig=1,ngrid
2050	if(tau_col(ig).gt.1.e3)then
2051	print*,'WARNING: tau_col=',tau_col(ig)
2052	print*,'at ig=',ig,'in PHYSIQ'
2053	endif
2054	end do
2055
2056	call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2057
2058	enddo
2059	endif
2060
2061	! output spectrum
2062	if(specOLR.and.corrk)then
2063	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2064	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2065	endif
2066
2067
2068	icount=icount+1
2069
2070	!!! DEALLOCATE STUFF
2071	if (lastcall) then
2072	IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom ) !! this was allocated in su_gases.F90
2073	IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac ) !! this was allocated in su_gases.F90
2074	endif
2075
2076
2077	return
2078	end subroutine physiq

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

Download in other formats:

Original Format