Context Navigation

physiq.F90 @ 374

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

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

Download in other formats:

Original Format