Context Navigation

physiq.F90 @ 304

Last change on this file since 304 was 253, checked in by emillour, 13 years ago

Generic GCM

Massive update to version 0.7

EM+RW

Property svn:executable set to *

File size: 77.4 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
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 pressure broadening broadening study
395	integer gascut
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 ice_min(ngridmx)
407	save ice_min
408	save ice_initial
409
410	integer num_run
411	logical ice_update
412	save ice_update
413
414	!=======================================================================
415
416
417	! 1. Initialisation
418	! -----------------
419
420	! 1.1 Initialisation only at first call
421	! ---------------------------------------
422	if (firstcall) then
423
424	if(ngrid.eq.1)then
425	saveG1D=day_step*10
426	!saveG1D=1
427	!saveG1D=-1 ! to never save
428	countG1D=1
429	endif
430
431	! variables set to 0
432	! ~~~~~~~~~~~~~~~~~~
433	dtrad(:,:) = 0.0
434	fluxrad(:) = 0.0
435	tau_col(:) = 0.0
436	zdtsw(:,:) = 0.0
437	zdtlw(:,:) = 0.0
438
439	! initialize tracer names, indexes and properties
440	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
441	tracerdyn=tracer
442	if (tracer) then
443	call initracer()
444	endif ! end tracer
445
446	! read startfi (initial parameters)
447	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
448
449	call phyetat0("startfi.nc",0,0,nsoilmx,nq, &
450	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
451	cloudfrac,totcloudfrac,hice)
452
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,6).eq.0)then
505	print*,'Updating ice at end of this year!'
506	ice_update=.true.
507	ice_min(:)=1.e4
508	endif
509	endif
510
511	! define surface as continent or ocean
512	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
513	! print*,'OCEAN EVERYWHERE FOR GJ581d RUNS'
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	watertest=.true. ! standard for now
567	if(ngrid.eq.1)then
568	qzero1D = 0.0
569	qsurf(1,igcm_h2o_vap) = qzero1D
570	do l=1, nlayer
571	pq(1,l,igcm_h2o_vap)=0.0
572	pq(1,l,igcm_h2o_ice)=0.0
573	enddo
574	endif
575
576	do ig=1,ngrid
577	qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
578	if(ice_update)then
579	ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
580	endif
581	enddo
582
583	endif
584	call su_watercycle ! even if we don't have a water cycle, we might
585	! need epsi for the wvp definitions in callcorrk.F
586
587	endif ! (end of "if firstcall")
588
589	! ---------------------------------------------------
590	! 1.2 Initializations done at every physical timestep:
591	! ---------------------------------------------------
592
593	if (ngrid.NE.ngridmx) then
594	print*,'STOP in PHYSIQ'
595	print*,'Probleme de dimensions :'
596	print*,'ngrid = ',ngrid
597	print*,'ngridmx = ',ngridmx
598	stop
599	endif
600
601	! Initialize various variables
602	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
603
604	pdu(:,:) = 0.0
605	pdv(:,:) = 0.0
606	! if ( (.not.nearco2cond).and.(.not.firstcall) ) then
607	if ( .not.nearco2cond ) then
608	pdt(:,:) = 0.0
609	endif ! this was the source of an evil bug...
610	pdq(:,:,:) = 0.0
611	pdpsrf(:) = 0.0
612	zflubid(:) = 0.0
613	zdtsurf(:) = 0.0
614	dqsurf(:,:) = 0.0
615
616	zday=pday+ptime ! compute time, in sols (and fraction thereof)
617
618	! Compute Stellar Longitude (Ls)
619	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
620	if (season) then
621	call stellarlong(zday,zls)
622	else
623	call stellarlong(float(day_ini),zls)
624	end if
625
626	! Compute geopotential between layers
627	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
628
629	do l=1,nlayer
630	do ig=1,ngrid
631	zzlay(ig,l)=pphi(ig,l)/g
632	enddo
633	enddo
634	do ig=1,ngrid
635	zzlev(ig,1)=0.
636	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
637	enddo
638	do l=2,nlayer
639	do ig=1,ngrid
640	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
641	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
642	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
643	enddo
644	enddo
645	! Potential temperature calculation may not be the same in physiq and dynamic...
646
647	! Compute potential temperature
648	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
649
650	do l=1,nlayer
651	do ig=1,ngrid
652	call calc_cpp3d(cppNI,rcpNI,pt(ig,l),pplay(ig,l))
653	cpp3D(ig,l) = cppNI
654	rcp3D(ig,l) = rcpNI
655	enddo
656	enddo
657
658	do l=1,nlayer
659	do ig=1,ngrid
660	zh(ig,l) = pt(ig,l)
661	zpopskNI(ig,l) = pt(ig,l)/zh(ig,l)
662	zpopsk(ig,l) = pt(ig,l)/zh(ig,l)
663	! we're only after zpopskNI here, not zh
664	! zh is calculated seperately before both vdifc and convadj
665	enddo
666	enddo
667
668	do l=1,nlayer
669	do ig=1,ngrid
670	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
671	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
672	enddo
673	enddo
674
675	!-----------------------------------------------------------------------
676	! 2. Compute radiative tendencies
677	!-----------------------------------------------------------------------
678
679	if (callrad) then
680	if( mod(icount-1,iradia).eq.0) then
681
682	! Compute local stellar zenith angles
683	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
684	call orbite(zls,dist_star,declin)
685
686	if (tlocked) then
687	ztim1=SIN(declin)
688	ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
689	ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
690
691	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
692	ztim1,ztim2,ztim3,mu0,fract)
693
694	elseif (diurnal) then
695	ztim1=SIN(declin)
696	ztim2=COS(declin)COS(2.pi*(zday-.5))
697	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
698
699	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
700	ztim1,ztim2,ztim3,mu0,fract)
701
702	else
703
704	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
705	! WARNING: this function appears not to work in 1D
706
707	endif
708
709	if (corrk) then
710
711	! a) Call correlated-k radiative transfer scheme
712	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
713
714	if(kastprof)then
715	ISR=fract(1)Fat1AU/dist_star*2/2.0
716
717	call simpleprof_fn(tsurf,pt,zpopsk,ISR,pplev(1,1))
718	open(116,file='profpres.out',form='formatted')
719	open(117,file='proftemp.out',form='formatted')
720	do l=1,nlayer
721	write(116,*) pplay(1,l)
722	write(117,*) pt(1,l)
723	enddo
724	close(116)
725	close(117)
726	endif
727
728
729	! Option to call scheme once for clear regions, once for cloudy regions (BC)
730	if(CLFvarying)then
731
732	clearsky=.true.
733	call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf, &
734	albedo,emis,mu0,pplev,pplay,pt, &
735	tsurf,fract,dist_star,aerosol,cpp3D, &
736	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
737	fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,ptime,pday, &
738	cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
739
740	clearsky=.false.
741	call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf, &
742	albedo,emis,mu0,pplev,pplay,pt, &
743	tsurf,fract,dist_star,aerosol,cpp3D, &
744	zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, &
745	fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,ptime,pday, &
746	cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
747
748	! Sum the fluxes and heating rates from cloudy/clear cases
749	do ig=1,ngrid
750	tf=totcloudfrac(ig)
751	ntf=1.-tf
752
753	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw2(ig)
754	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw2(ig)
755	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw2(ig)
756	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw2(ig)
757	tau_col(ig) = ntftau_col1(ig) + tftau_col2(ig)
758
759	do l=1,nlayer
760	zdtlw(ig,l) = ntfzdtlw1(ig,l) + tfzdtlw2(ig,l)
761	zdtsw(ig,l) = ntfzdtsw1(ig,l) + tfzdtsw2(ig,l)
762	enddo
763	enddo
764
765	! Otherwise standard callcorrk
766	else
767
768	clearsky=.false.
769	call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf, &
770	albedo,emis,mu0,pplev,pplay,pt, &
771	tsurf,fract,dist_star,aerosol,cpp3D, &
772	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
773	fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
774	cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
775
776	endif
777
778	! Radiative flux from the sky absorbed by the surface (W.m-2)
779	GSR=0.0
780	do ig=1,ngrid
781	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig) &
782	+fluxsurf_sw(ig)*(1.-albedo(ig))
783
784	if(noradsurf)then ! no lower surface; SW flux just disappears
785	GSR = GSR + fluxsurf_sw(ig)*area(ig)
786	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
787	endif
788
789	enddo
790	if(noradsurf)then
791	print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
792	endif
793
794	! Net atmospheric radiative heating rate (K.s-1)
795	do l=1,nlayer
796	do ig=1,ngrid
797	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
798	enddo
799	enddo
800
801
802	elseif(newtonian)then
803
804	! b) Call Newtonian cooling scheme
805	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806	call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
807
808	do ig=1,ngrid
809	zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
810	enddo
811	! e.g. surface becomes proxy for 1st atmospheric layer ?
812
813	else
814
815	! c) Atmosphere has no radiative effect
816	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
817	do ig=1,ngrid
818	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_star**2
819	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
820	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
821	endif
822	fluxsurf_sw(ig) = fluxtop_dn(ig)
823	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
824	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
825	enddo ! radiation skips the atmosphere entirely
826
827	do l=1,nlayer
828	do ig=1,ngrid
829	dtrad(ig,l)=0.0
830	enddo
831	enddo ! hence no atmospheric radiative heating
832
833	endif ! if corrk
834
835	endif ! of if(mod(icount-1,iradia).eq.0)
836
837
838	! Transformation of the radiative tendencies
839	! ------------------------------------------
840
841	do ig=1,ngrid
842	zplanck(ig)=tsurf(ig)*tsurf(ig)
843	zplanck(ig)=emis(ig)stephanzplanck(ig)*zplanck(ig)
844	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
845	enddo
846
847	do l=1,nlayer
848	do ig=1,ngrid
849	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
850	enddo
851	enddo
852
853	!-------------------------
854	! test energy conservation
855	if(enertest)then
856	dEtotSW = 0.0
857	dEtotLW = 0.0
858	dEtotsSW = 0.0
859	dEtotsLW = 0.0
860	do ig = 1, ngrid
861	do l = 1, nlayer
862	masse = (pplev(ig,l) - pplev(ig,l+1))/g
863	dEtotSW = dEtotSW + cpp3D(ig,l)massezdtsw(ig,l)*area(ig)
864	dEtotLW = dEtotLW + cpp3D(ig,l)massezdtlw(ig,l)*area(ig)
865	enddo
866	dEtotsSW = dEtotsSW + fluxsurf_sw(ig)(1.-albedo(ig))area(ig)
867	dEtotsLW = dEtotsLW + emis(ig)fluxsurf_lw(ig)area(ig)-zplanck(ig)*area(ig)
868	enddo
869	dEtotSW = dEtotSW/totarea
870	dEtotLW = dEtotLW/totarea
871	dEtotsSW = dEtotsSW/totarea
872	dEtotsLW = dEtotsLW/totarea
873	print*,'---------------------------------------------------------------'
874	print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2'
875	print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2'
876	print*,'In corrk SW surface energy change =',dEtotsSW,' W m-2'
877	print*,'In corrk LW surface energy change =',dEtotsLW,' W m-2'
878	endif
879	!-------------------------
880
881	endif ! of if (callrad)
882
883
884	!-----------------------------------------------------------------------
885	! 4. Vertical diffusion (turbulent mixing):
886	! -----------------------------------------
887
888	if (calldifv) then
889
890	do ig=1,ngrid
891	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
892	enddo
893
894	zdum1(:,:)=0.0
895	zdum2(:,:)=0.0
896	do l=1,nlayer
897	do ig=1,ngrid
898	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
899	enddo
900	enddo
901
902	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
903	ptimestep,capcal,lwrite, &
904	pplay,pplev,zzlay,zzlev,z0, &
905	pu,pv,zh,pq,tsurf,emis,qsurf, &
906	zdum1,zdum2,zdh,pdq,zflubid, &
907	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
908	zdqdif,zdqsdif)
909
910	do l=1,nlayer
911	do ig=1,ngrid
912	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
913	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
914	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
915	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
916	enddo
917	enddo
918
919	do ig=1,ngrid
920	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
921	enddo
922
923	if (tracer) then
924	do iq=1, nq
925	do l=1,nlayer
926	do ig=1,ngrid
927	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
928	enddo
929	enddo
930	enddo
931	do iq=1, nq
932	do ig=1,ngrid
933	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
934	enddo
935	enddo
936
937	end if ! of if (tracer)
938
939	!-------------------------
940	! test energy conservation
941	if(enertest)then
942	dEtot=0.0
943	dEtots=0.0
944	vdifcdE(:)=0.0
945	do ig = 1, ngrid
946	do l = 1, nlayer
947	masse = (pplev(ig,l) - pplev(ig,l+1))/g
948	dEtot = dEtot + cpp3D(ig,l)massezdtdif(ig,l)*area(ig)
949
950	vabs = sqrt(pdu(ig,l)2 + pdv(ig,l)2)
951	dvabs = sqrt(zdudif(ig,l)2 + zdvdif(ig,l)2)
952	dEtot = dEtot + massevabsdvabs*area(ig)
953
954	vdifcdE(ig) = vdifcdE(ig) + massevabsdvabs
955
956	enddo
957	dEtots = dEtots + capcal(ig)zdtsdif(ig)area(ig)
958	vdifcdE(ig) = vdifcdE(ig) + capcal(ig)*zdtsdif(ig)
959	enddo
960	dEtot = dEtot/totarea
961	dEtots = dEtots/totarea
962	print*,'In difv atmospheric energy change =',dEtot,' W m-2'
963	print*,'In difv surface energy change =',dEtots,' W m-2'
964	!print*,'Note we ignore the wind change...'
965	! and latent heat for that matter...
966	endif
967	!-------------------------
968
969	!-------------------------
970	! test water conservation
971	if(watertest.and.water)then
972	dWtot=0.0
973	dWtots=0.0
974	nconsMAX=0.0
975	do ig = 1, ngrid
976
977	vdifcncons(ig)=0.0
978	do l = 1, nlayer
979	masse = (pplev(ig,l) - pplev(ig,l+1))/g
980
981	iq = igcm_h2o_vap
982	dWtot = dWtot + massezdqdif(ig,l,iq)ptimestep*area(ig)
983	vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
984
985	iq = igcm_h2o_ice
986	dWtot = dWtot + massezdqdif(ig,l,iq)ptimestep*area(ig)
987	vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
988
989	enddo
990
991	iq = igcm_h2o_vap
992	dWtots = dWtots + zdqsdif(ig,iq)ptimesteparea(ig)
993	vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
994
995	iq = igcm_h2o_ice
996	dWtots = dWtots + zdqsdif(ig,iq)ptimesteparea(ig)
997	vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
998
999	if(vdifcncons(ig).gt.nconsMAX)then
1000	nconsMAX=vdifcncons(ig)
1001	endif
1002
1003	enddo
1004
1005	dWtot = dWtot/totarea
1006	dWtots = dWtots/totarea
1007	print*,'---------------------------------------------------------------'
1008	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
1009	print*,'In difv surface water change =',dWtots,' kg m-2'
1010	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
1011	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
1012
1013
1014	endif
1015	!-------------------------
1016
1017	else
1018
1019	if(.not.newtonian)then
1020
1021	do ig=1,ngrid
1022	zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
1023	enddo
1024
1025	endif
1026
1027	endif ! of if (calldifv)
1028
1029
1030	!-----------------------------------------------------------------------
1031	! 5. Dry convective adjustment:
1032	! -----------------------------
1033
1034	if(calladj) then
1035
1036	do l=1,nlayer
1037	do ig=1,ngrid
1038	if(nonideal)then
1039	zdh(ig,l) = pdt(ig,l)/zpopskNI(ig,l)
1040	else
1041	zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
1042	endif
1043	enddo
1044	enddo
1045	zduadj(:,:)=0.0
1046	zdvadj(:,:)=0.0
1047	zdhadj(:,:)=0.0
1048
1049	if(nonideal)then
1050	print*,'Nonideal is not used at the moment'
1051	call abort
1052	call convadj(ngrid,nlayer,nq,ptimestep, &
1053	pplay,pplev,zpopskNI, &
1054	pu,pv,zh,pq, &
1055	pdu,pdv,zdh,pdq, &
1056	zduadj,zdvadj,zdhadj, &
1057	zdqadj)
1058	else
1059
1060	call convadj(ngrid,nlayer,nq,ptimestep, &
1061	pplay,pplev,zpopsk, &
1062	pu,pv,zh,pq, &
1063	pdu,pdv,zdh,pdq, &
1064	zduadj,zdvadj,zdhadj, &
1065	zdqadj)
1066
1067	endif
1068
1069	do l=1,nlayer
1070	do ig=1,ngrid
1071	pdu(ig,l) = pdu(ig,l) + zduadj(ig,l)
1072	pdv(ig,l) = pdv(ig,l) + zdvadj(ig,l)
1073	if(nonideal)then
1074	pdt(ig,l) = pdt(ig,l) + zdhadj(ig,l)*zpopskNI(ig,l)
1075	zdtadj(ig,l) = zdhadj(ig,l)*zpopskNI(ig,l) ! for diagnostic only
1076	else
1077	pdt(ig,l) = pdt(ig,l) + zdhadj(ig,l)*zpopsk(ig,l)
1078	zdtadj(ig,l) = zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1079	endif
1080	enddo
1081	enddo
1082
1083	if(tracer) then
1084	do iq=1, nq
1085	do l=1,nlayer
1086	do ig=1,ngrid
1087	pdq(ig,l,iq) = pdq(ig,l,iq) + zdqadj(ig,l,iq)
1088	enddo
1089	enddo
1090	enddo
1091	end if
1092
1093	!-------------------------
1094	! test energy conservation
1095	if(enertest)then
1096	dEtot=0.0
1097	do ig = 1, ngrid
1098	do l = 1, nlayer
1099	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1100	dEtot = dEtot + cpp3D(ig,l)massezdtadj(ig,l)*area(ig)
1101	enddo
1102	enddo
1103	dEtot=dEtot/totarea
1104	print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1105	endif
1106	!-------------------------
1107
1108	!-------------------------
1109	! test water conservation
1110	if(watertest)then
1111	dWtot=0.0
1112	do ig = 1, ngrid
1113	cadjncons(ig)=0.0
1114	do l = 1, nlayer
1115	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1116
1117	iq = igcm_h2o_vap
1118	dWtot = dWtot + massezdqadj(ig,l,iq)ptimestep*area(ig)
1119	cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
1120
1121	iq = igcm_h2o_ice
1122	dWtot = dWtot + massezdqadj(ig,l,iq)ptimestep*area(ig)
1123	cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
1124
1125	enddo
1126
1127	enddo
1128	dWtot=dWtot/totarea
1129	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
1130	endif
1131	!-------------------------
1132
1133	endif ! of if(calladj)
1134
1135	!-----------------------------------------------------------------------
1136	! 6. Carbon dioxide condensation-sublimation:
1137	! -------------------------------------------
1138
1139	if (co2cond) then
1140	if (.not.tracer) then
1141	print*,'We need a CO2 ice tracer to condense CO2'
1142	call abort
1143	endif
1144
1145	call condens_co2cloud(ngrid,nlayer,nq,ptimestep, &
1146	capcal,pplay,pplev,tsurf,pt, &
1147	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1148	qsurf(1,igcm_co2_ice),albedo,emis, &
1149	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1150	zdqc,reffrad,cpp3D)
1151
1152	do l=1,nlayer
1153	do ig=1,ngrid
1154	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1155	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1156	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1157	enddo
1158	enddo
1159	do ig=1,ngrid
1160	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1161	enddo
1162
1163	do iq=1,nq ! should use new notation here !
1164	do l=1,nlayer
1165	do ig=1,ngrid
1166	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1167	enddo
1168	enddo
1169	enddo
1170	! Note: we do not add surface co2ice tendency
1171	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
1172
1173	!-------------------------
1174	! test energy conservation
1175	if(enertest)then
1176	dEtot=0.0
1177	dEtots=0.0
1178	do ig = 1, ngrid
1179	do l = 1, nlayer
1180	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1181	dEtot = dEtot + cpp3D(ig,l)massezdtc(ig,l)*area(ig)
1182	enddo
1183	dEtots = dEtots + capcal(ig)zdtsurfc(ig)area(ig)
1184	enddo
1185	dEtot=dEtot/totarea
1186	dEtots=dEtots/totarea
1187	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
1188	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
1189	endif
1190	!-------------------------
1191
1192	endif ! of if (co2cond)
1193
1194
1195	!-----------------------------------------------------------------------
1196	! 7. Specific parameterizations for tracers
1197	! -----------------------------------------
1198
1199	if (tracer) then
1200
1201	! 7a. Water and ice
1202	! ---------------
1203	if (water) then
1204
1205	! ----------------------------------------
1206	! Water ice condensation in the atmosphere
1207	! ----------------------------------------
1208	if(watercond)then
1209
1210	if(RLVTT.gt.1.e-8)then
1211
1212	! Re-evaporate cloud water/ice
1213	call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1214	DO l = 1, nlayer
1215	DO ig = 1, ngrid
1216	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1217	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1218	pdt(ig,l) = pdt(ig,l)+dtevap(ig,l)
1219	enddo
1220	enddo
1221
1222	call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1223	do l=1,nlayer
1224	do ig=1,ngrid
1225	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
1226	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
1227	pdt(ig,l) = pdt(ig,l)+dtmoist(ig,l)
1228	enddo
1229	enddo
1230
1231	!-------------------------
1232	! test energy conservation
1233	if(enertest)then
1234	dEtot=0.0
1235	madjdE(:)=0.0
1236	do ig = 1, ngrid
1237	do l = 1, nlayer
1238	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1239	dEtot = dEtot + cpp3D(ig,l)masse(dtmoist(ig,l)+dtevap(ig,l))*area(ig)
1240	madjdE(ig) = madjdE(ig) + cpp3D(ig,l)masse(dtmoist(ig,l)+dtevap(ig,l))
1241	enddo
1242	enddo
1243	dEtot=dEtot/totarea
1244	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
1245	endif
1246	!-------------------------
1247
1248	!-------------------------
1249	! test water conservation
1250	if(watertest)then
1251	dWtot=0.0
1252	do ig = 1, ngrid
1253	do iq = 1 , nq
1254	do l = 1, nlayer
1255	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1256	dWtot = dWtot + massedqmoist(ig,l,igcm_h2o_vap)area(ig)*ptimestep
1257	dWtot = dWtot + massedqmoist(ig,l,igcm_h2o_ice)area(ig)*ptimestep
1258	enddo
1259	enddo
1260	enddo
1261	dWtot=dWtot/totarea
1262	print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
1263	endif
1264	!-------------------------
1265
1266
1267	endif
1268
1269
1270	! Re-evaporate cloud water/ice
1271	call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1272	do l = 1, nlayer
1273	do ig = 1, ngrid
1274	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1275	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1276	pdt(ig,l) = pdt(ig,l)+dtevap(ig,l)
1277	enddo
1278	enddo ! note: we use qevap but not tevap in largescale/moistadj
1279	! otherwise is a big mess
1280
1281	call largescale(ptimestep,pplev,pplay,pt,qevap, & ! a bug was here!
1282	pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
1283	do l=1,nlayer
1284	do ig=1,ngrid
1285	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
1286	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
1287	pdt(ig,l) = pdt(ig,l)+dtlscale(ig,l)
1288
1289	if(.not.aerofixed)then
1290	reffrad(ig,l,2)=reffH2O(ig,l)
1291	endif
1292
1293	enddo
1294	enddo
1295
1296	!-------------------------
1297	! test energy conservation
1298	if(enertest)then
1299	dEtot=0.0
1300	lscaledE(:)=0.0
1301	do ig = 1, ngrid
1302	do l = 1, nlayer
1303	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1304	dEtot = dEtot + cpp3D(ig,l)masse(dtlscale(ig,l)+dtevap(ig,l))*area(ig)
1305	lscaledE(ig) = lscaledE(ig) + cpp3D(ig,l)masse(dtlscale(ig,l)+dtevap(ig,l))
1306	enddo
1307	enddo
1308	dEtot=dEtot/totarea
1309	print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
1310	endif
1311	!-------------------------
1312
1313	!-------------------------
1314	! test water conservation
1315	if(watertest)then
1316	dWtot=0.0
1317	do ig = 1, ngrid
1318	do iq = 1 , nq
1319	do l = 1, nlayer
1320	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1321	dWtot = dWtot + massedqvaplscale(ig,l)area(ig)*ptimestep
1322	dWtot = dWtot + massedqcldlscale(ig,l)area(ig)*ptimestep
1323	enddo
1324	enddo
1325	enddo
1326	dWtot=dWtot/totarea
1327	print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
1328	endif
1329	!-------------------------
1330
1331	! compute cloud fraction
1332	do l = 1, nlayer
1333	do ig = 1,ngrid
1334	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1335	enddo
1336	enddo
1337
1338	! compute total cloud fraction in column
1339	call totalcloudfrac(cloudfrac,totcloudfrac)
1340
1341	endif ! of if (watercondense)
1342
1343
1344	! --------------------------------
1345	! Water ice / liquid precipitation
1346	! --------------------------------
1347	if(waterrain)then
1348
1349	zdqrain(:,:,:) = 0.0
1350	zdqsrain(:) = 0.0
1351	zdqssnow(:) = 0.0
1352
1353	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
1354	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
1355
1356	do l=1,nlayer
1357	do ig=1,ngrid
1358	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
1359	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
1360	pdt(ig,l) = pdt(ig,l)+zdtrain(ig,l)
1361	enddo
1362	enddo
1363
1364	do ig=1,ngrid
1365	dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
1366	dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
1367	rainout(ig) = zdqsrain(ig)+zdqssnow(ig) ! diagnostic
1368	enddo
1369
1370	!-------------------------
1371	! test energy conservation
1372	if(enertest)then
1373	dEtot=0.0
1374	do ig = 1, ngrid
1375	do l = 1, nlayer
1376	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1377	dEtot = dEtot + cpp3D(ig,l)massezdtrain(ig,l)*area(ig)
1378	enddo
1379	enddo
1380	dEtot=dEtot/totarea
1381	print*,'In rain atmospheric energy change =',dEtot,' W m-2'
1382	endif
1383	!-------------------------
1384
1385
1386	!-------------------------
1387	! test energy conservation
1388	if(enertest)then
1389	dEtot=0.0
1390	do ig = 1, ngrid
1391	do l = 1, nlayer
1392	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1393	dEtot = dEtot + cpp3D(ig,l)massezdtrain(ig,l)*area(ig)
1394	enddo
1395	enddo
1396	dEtot=dEtot/totarea
1397	print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
1398
1399	dEtot=0.0
1400	do ig = 1, ngrid
1401	do l = 1, nlayer
1402	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1403	dItot = dItot + massezdqrain(ig,l,igcm_h2o_ice)area(ig)
1404	dVtot = dVtot + massezdqrain(ig,l,igcm_h2o_vap)area(ig)
1405	enddo
1406	dItot = dItot + zdqssnow(ig)*area(ig)
1407	dVtot = dVtot + zdqsrain(ig)*area(ig)
1408	enddo
1409	dEtot=(dItotRLVTT/cpp + dVtotRLVTT/cpp)/totarea
1410	print,'In rain dItot =',dItotRLVTT/(cpp*totarea),' W m-2'
1411	print,'In rain dVtot =',dVtotRLVTT/(cpp*totarea),' W m-2'
1412	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
1413	endif
1414	!-------------------------
1415
1416	!-------------------------
1417	! test water conservation
1418	if(watertest)then
1419	dWtot=0.0
1420	dWtots=0.0
1421	do ig = 1, ngrid
1422	!do iq = 1 , nq
1423	do l = 1, nlayer
1424	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1425	dWtot = dWtot + massezdqrain(ig,l,igcm_h2o_vap)area(ig)*ptimestep
1426	dWtot = dWtot + massezdqrain(ig,l,igcm_h2o_ice)area(ig)*ptimestep
1427	enddo
1428	!enddo
1429	dWtots = dWtots + (zdqsrain(ig)+zdqssnow(ig))area(ig)ptimestep
1430	enddo
1431	dWtot=dWtot/totarea
1432	dWtots=dWtots/totarea
1433	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
1434	print*,'In rain surface water change =',dWtots,' kg m-2'
1435	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
1436	endif
1437	!-------------------------
1438
1439	end if ! of if (waterrain)
1440	end if ! of if (water)
1441
1442
1443	! 7c. Aerosol particles
1444	! -------------------
1445	! -------------
1446	! Sedimentation
1447	! -------------
1448	if (sedimentation) then
1449	zdqsed(:,:,:) = 0.0
1450	zdqssed(:,:) = 0.0
1451
1452
1453	!-------------------------
1454	! find qtot
1455	if(watertest)then
1456	dWtot=0.0
1457	dWtots=0.0
1458	iq=3
1459	do ig = 1, ngrid
1460	do l = 1, nlayer
1461	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1462	dWtot = dWtot + massepq(ig,l,iq)area(ig)*ptimestep
1463	dWtots = dWtots + massepdq(ig,l,iq)area(ig)*ptimestep
1464	enddo
1465	enddo
1466	dWtot=dWtot/totarea
1467	dWtots=dWtots/totarea
1468	print*,'Before sedim pq =',dWtot,' kg m-2'
1469	print*,'Before sedim pdq =',dWtots,' kg m-2'
1470	endif
1471	!-------------------------
1472
1473	call callsedim(ngrid,nlayer,ptimestep, &
1474	pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
1475
1476	!-------------------------
1477	! find qtot
1478	if(watertest)then
1479	dWtot=0.0
1480	dWtots=0.0
1481	iq=3
1482	do ig = 1, ngrid
1483	do l = 1, nlayer
1484	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1485	dWtot = dWtot + massepq(ig,l,iq)area(ig)*ptimestep
1486	dWtots = dWtots + massepdq(ig,l,iq)area(ig)*ptimestep
1487	enddo
1488	enddo
1489	dWtot=dWtot/totarea
1490	dWtots=dWtots/totarea
1491	print*,'After sedim pq =',dWtot,' kg m-2'
1492	print*,'After sedim pdq =',dWtots,' kg m-2'
1493	endif
1494	!-------------------------
1495
1496	do iq=1,nq
1497	! for now, we only allow H2O ice to sediment
1498	! and as in rain.F90, whether it falls as rain or snow depends
1499	! only on the surface temperature
1500	do ig=1,ngrid
1501	do l=1,nlayer
1502	pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
1503	enddo
1504	dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
1505	enddo
1506	enddo
1507
1508	!-------------------------
1509	! test water conservation
1510	if(watertest)then
1511	dWtot=0.0
1512	dWtots=0.0
1513	do iq=1,nq
1514	do ig = 1, ngrid
1515	do l = 1, nlayer
1516	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1517	dWtot = dWtot + massezdqsed(ig,l,iq)area(ig)*ptimestep
1518	enddo
1519	dWtots = dWtots + zdqssed(ig,iq)area(ig)ptimestep
1520	enddo
1521	enddo
1522	dWtot=dWtot/totarea
1523	dWtots=dWtots/totarea
1524	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
1525	print*,'In sedim surface ice change =',dWtots,' kg m-2'
1526	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
1527	endif
1528	!-------------------------
1529
1530	end if ! of if (sedimentation)
1531
1532
1533	! 7d. Updates
1534	! ---------
1535
1536	! ---------------------------------
1537	! Updating tracer budget on surface
1538	! ---------------------------------
1539
1540	if(hydrology)then
1541
1542	call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
1543	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
1544	! note: for now, also changes albedo in the subroutine
1545
1546	do ig=1,ngrid
1547	zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
1548	do iq=1,nq
1549	qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
1550	enddo
1551	enddo
1552	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
1553
1554	!-------------------------
1555	! test energy conservation
1556	if(enertest)then
1557	dEtot=0.0
1558	do ig = 1, ngrid
1559	dEtots = dEtots + capcal(ig)zdtsurf_hyd(ig)area(ig)
1560	enddo
1561	dEtot=dEtot/totarea
1562	print*,'In hydrol atmospheric energy change =',dEtot,' W m-2'
1563	endif
1564	!-------------------------
1565
1566	!-------------------------
1567	! test water conservation
1568	if(watertest)then
1569	dWtots=0.0
1570	do ig = 1, ngrid
1571	dWtots = dWtots + dqs_hyd(ig,igcm_h2o_ice)area(ig)ptimestep
1572	enddo
1573	dWtots=dWtots/totarea
1574	print*,'In hydrol surface ice change =',dWtots,' kg m-2'
1575	dWtots=0.0
1576	do ig = 1, ngrid
1577	dWtots = dWtots + dqs_hyd(ig,igcm_h2o_vap)area(ig)ptimestep
1578	enddo
1579	dWtots=dWtots/totarea
1580	print*,'In hydrol surface water change =',dWtots,' kg m-2'
1581	print*,'---------------------------------------------------------------'
1582	endif
1583	!-------------------------
1584
1585	ELSE ! of if (hydrology)
1586
1587	do iq=1,nq
1588	do ig=1,ngrid
1589	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1590	enddo
1591	enddo
1592
1593	END IF ! of if (hydrology)
1594
1595	! Add qsurf to qsurf_hist, which is what we save in
1596	! diagfi.nc etc. At the same time, we set the water
1597	! content of ocean gridpoints back to zero, in order
1598	! to avoid rounding errors in vdifc, rain
1599	do ig = 1, ngrid
1600	do iq = 1, nq
1601	if(iq.eq.igcm_h2o_vap .and. rnat(ig).eq.0)then ! if liquid water and terrain = ocean
1602	qsurf_hist(ig,iq) = qsurf(ig,iq)
1603	!qsurf(ig,iq) = qcol(ig,iq)
1604	! the value of qsurf we choose here makes NO DIFFERENCE TO ANYTHING AT ALL
1605	else
1606	qsurf_hist(ig,iq) = qsurf(ig,iq)
1607	endif
1608	enddo
1609	enddo
1610
1611	if(ice_update)then
1612	do ig = 1, ngrid
1613	ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
1614	enddo
1615	endif
1616
1617	endif ! of if (tracer)
1618
1619
1620	! print*,'we got this far 1'
1621
1622	!-----------------------------------------------------------------------
1623	! 9. Surface and sub-surface soil temperature
1624	!-----------------------------------------------------------------------
1625
1626
1627	! Increment surface temperature
1628	if(.not.kastprof)then
1629	do ig=1,ngrid
1630	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1631	enddo
1632	endif
1633
1634	! Compute soil temperatures and subsurface heat flux
1635	if (callsoil) then
1636	call soil(ngrid,nsoilmx,.false.,inertiedat, &
1637	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1638	endif
1639
1640	!-------------------------
1641	! test energy conservation
1642	if(enertest)then
1643	dEtots=0.0
1644	do ig = 1, ngrid
1645	dEtots = dEtots + capcal(ig)zdtsurf(ig)area(ig)
1646	enddo
1647	dEtots=dEtots/totarea
1648	print*,'Surface energy change=',dEtots,' W m-2'
1649	endif
1650	!-------------------------
1651
1652	!-----------------------------------------------------------------------
1653	! 10. Perform diagnostics and write output files
1654	!-----------------------------------------------------------------------
1655
1656	! -------------------------------
1657	! Dynamical fields incrementation
1658	! -------------------------------
1659	! For output only: the actual model integration is performed in the dynamics
1660
1661	! temperature, zonal and meridional wind
1662	do l=1,nlayer
1663	do ig=1,ngrid
1664	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
1665	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
1666	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
1667
1668	! diagnostic
1669	zdtdyn(ig,l) = ztprevious(ig,l)-pt(ig,l)
1670	ztprevious(ig,l) = zt(ig,l)
1671	enddo
1672	enddo
1673
1674	if(firstcall)then
1675	zdtdyn(:,:)=0.0
1676	endif
1677
1678	! dynamical heating diagnostic
1679	fluxdyn(:)=0.
1680	do ig=1,ngrid
1681	do l=1,nlayer
1682	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep) &
1683	(pplev(ig,l)-pplev(ig,l+1))cpp3D(ig,l)/g
1684	enddo
1685	enddo
1686
1687	! tracers
1688	do iq=1, nq
1689	do l=1,nlayer
1690	do ig=1,ngrid
1691	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
1692	enddo
1693	enddo
1694	enddo
1695
1696	! surface pressure
1697	do ig=1,ngrid
1698	ps(ig) = pplev(ig,1) + pdpsrf(ig)*ptimestep
1699	enddo
1700
1701	! pressure
1702	do l=1,nlayer
1703	do ig=1,ngrid
1704	zplev(ig,l) = pplev(ig,l)/pplev(ig,1)*ps(ig)
1705	zplay(ig,l) = pplay(ig,l)/pplev(ig,1)*ps(ig)
1706	enddo
1707	enddo
1708
1709	! ---------------------------------------------------------
1710	! Surface and soil temperature information
1711	! ---------------------------------------------------------
1712
1713	Ts1 = 0.0
1714	Ts2 = 99999.9
1715	Ts3 = 0.0
1716	TsS = 0.0 ! mean temperature at bottom soil layer
1717	do ig=1,ngrid
1718	Ts1 = Ts1 + area(ig)*tsurf(ig)
1719	Ts2 = min(Ts2,tsurf(ig))
1720	Ts3 = max(Ts3,tsurf(ig))
1721	TsS = TsS + area(ig)*tsoil(ig,nsoilmx)
1722	end do
1723	Ts1=Ts1/totarea
1724	TsS=TsS/totarea
1725	if(callsoil)then
1726	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1727	print*,Ts1,Ts2,Ts3,TsS
1728	else
1729	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1730	print*,Ts1,Ts2,Ts3
1731	endif
1732
1733	! ---------------------------------------------------------
1734	! Check the energy balance of the simulation during the run
1735	! ---------------------------------------------------------
1736
1737	if(corrk)then
1738
1739	ISR = 0.0
1740	ASR = 0.0
1741	OLR = 0.0
1742	GND = 0.0
1743	DYN = 0.0
1744	do ig=1,ngrid
1745	ISR = ISR + area(ig)*fluxtop_dn(ig)
1746	ASR = ASR + area(ig)*fluxabs_sw(ig)
1747	OLR = OLR + area(ig)*fluxtop_lw(ig)
1748	GND = GND + area(ig)*fluxgrd(ig)
1749	DYN = DYN + area(ig)*fluxdyn(ig)
1750
1751	if(fluxtop_dn(ig).lt.0.0)then
1752	print*,'fluxtop_dn has gone crazy'
1753	print*,'fluxtop_dn=',fluxtop_dn(ig)
1754	print*,'tau_col=',tau_col(ig)
1755	print*,'aerosol=',aerosol(ig,:,:)
1756	print*,'temp= ',pt(ig,:)
1757	print*,'pplay= ',pplay(ig,:)
1758	call abort
1759	endif
1760	end do
1761
1762	if(ngridmx.eq.1)then
1763	DYN=0.0
1764	endif
1765
1766	print*,' ISR ASR OLR GND DYN [W m^-2]'
1767	print*, ISR/totarea,ASR/totarea,OLR/totarea,GND/totarea,DYN/totarea
1768
1769	if(enertest)then
1770	print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
1771	endif
1772
1773	if(meanOLR)then
1774	if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
1775	! to record global radiative balance
1776	open(92,file="rad_bal.out",form='formatted',access='append')
1777	write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
1778	close(92)
1779	open(93,file="tem_bal.out",form='formatted',access='append')
1780	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1781	close(93)
1782	endif
1783	endif
1784
1785	if(kastprof)then ! for f(p,T) computation (not needed in universal model)
1786	open(15,file='ene_bal.out',form='formatted',access='append')
1787	write(15,*) pplev(1,1), Tsurf(1), ISR, OLR, ASR
1788	close(15)
1789	open(16,file='Seff.out',form='formatted')
1790	write(16,*) OLR/ASR
1791	close(16)
1792	open(17,file='alb.out',form='formatted')
1793	write(17,*) 1.0-ASR/ISR
1794	close(17)
1795	call abort
1796	endif
1797
1798	endif
1799
1800	! ------------------------------------------------------------------
1801	! Diagnostic to test radiative-convective timescales in code
1802	! ------------------------------------------------------------------
1803	if(testradtimes)then
1804	open(38,file="tau_phys.out",form='formatted',access='append')
1805	ig=1
1806	do l=1,nlayer
1807	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1808	enddo
1809	close(38)
1810	print*,'As testradtimes enabled, exiting physics on first call'
1811	call abort
1812	endif
1813
1814	! ---------------------------------------------------------
1815	! Compute column amounts (kg m-2) if tracers are enabled
1816	! ---------------------------------------------------------
1817	if(tracer)then
1818	qcol(:,:)=0.0
1819	do iq=1,nq
1820	do ig=1,ngrid
1821	do l=1,nlayer
1822	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) * &
1823	(pplev(ig,l) - pplev(ig,l+1)) / g
1824	enddo
1825	enddo
1826	enddo
1827
1828	! not generalised for arbitrary aerosols yet!!!
1829	reffcol(:,:)=0.0
1830	do ig=1,ngrid
1831	do l=1,nlayer
1832	reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * reffrad(ig,l,1) * &
1833	(pplev(ig,l) - pplev(ig,l+1)) / g
1834	if(water)then
1835	reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * reffrad(ig,l,2) * &
1836	(pplev(ig,l) - pplev(ig,l+1)) / g
1837
1838	endif
1839	enddo
1840	enddo
1841
1842	endif
1843
1844	! ---------------------------------------------------------
1845	! Test for water conservation if water is enabled
1846	! ---------------------------------------------------------
1847
1848	if(water)then
1849
1850	icesrf = 0.0
1851	liqsrf = 0.0
1852	icecol = 0.0
1853	vapcol = 0.0
1854
1855	h2otot = 0.0
1856	do ig=1,ngrid
1857
1858	icesrf = icesrf + area(ig)*qsurf_hist(ig,igcm_h2o_ice)
1859	liqsrf = liqsrf + area(ig)*qsurf_hist(ig,igcm_h2o_vap)
1860	icecol = icecol + area(ig)*qcol(ig,igcm_h2o_ice)
1861	vapcol = vapcol + area(ig)*qcol(ig,igcm_h2o_vap)
1862
1863	h2otot = h2otot + area(ig)* &
1864	(qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap) &
1865	+qsurf_hist(ig,igcm_h2o_ice)+qsurf_hist(ig,igcm_h2o_vap))
1866	end do
1867
1868	print*,' Total water amount [kg m^-2]: ',h2otot/totarea
1869	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
1870	print*, icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
1871
1872	if(meanOLR)then
1873	if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
1874	! to record global water balance
1875	open(98,file="h2o_bal.out",form='formatted',access='append')
1876	write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
1877	close(98)
1878	endif
1879	endif
1880
1881	endif
1882
1883	! ---------------------------------------------------------
1884	! Calculate RH for diagnostic if water is enabled
1885	! ---------------------------------------------------------
1886
1887	if(water)then
1888	do l = 1, nlayer
1889	do ig = 1, ngrid
1890	call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
1891	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
1892	enddo
1893	enddo
1894
1895	! compute maximum possible H2O column amount (100% saturation)
1896	do ig=1,ngrid
1897	H2Omaxcol(ig)=0.0
1898	do l=1,nlayer
1899	H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) * &
1900	(pplev(ig,l) - pplev(ig,l+1))/g
1901	enddo
1902	enddo
1903
1904	endif
1905
1906
1907	if (ngrid.ne.1) then
1908	print*,''
1909	print,'--> Ls =',zls180./pi
1910	! -------------------------------------------------------------------
1911	! Writing NetCDF file "RESTARTFI" at the end of the run
1912	! -------------------------------------------------------------------
1913	! Note: 'restartfi' is stored just before dynamics are stored
1914	! in 'restart'. Between now and the writting of 'restart',
1915	! there will have been the itau=itau+1 instruction and
1916	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1917	! thus we store for time=time+dtvr
1918
1919	if(lastcall) then
1920	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1921
1922
1923	! Update surface ice distribution to iterate to steady state if requested
1924	if(ice_update)then
1925	do ig = 1, ngrid
1926
1927	! if ice amount is declining, set it to zero
1928	if(qsurf_hist(ig,igcm_h2o_ice).lt.ice_initial(ig))then
1929	qsurf_hist(ig,igcm_h2o_ice)=0.0
1930	qsurf_hist(ig,igcm_h2o_vap)=0.0
1931	endif
1932
1933	! if ice amount is increasing and was never less than 0.5 kg/m2,
1934	! set value in gridbox to 1x10^3 kg/m2
1935	!if(qsurf_hist(ig,igcm_h2o_ice).gt.0.0 .and. ice_min(ig).gt.0.5)then
1936	if(qsurf_hist(ig,igcm_h2o_ice).gt.ice_initial(ig))then
1937	if(ice_min(ig).gt.0.5)then ! ignore seasonal ice
1938	qsurf_hist(ig,igcm_h2o_ice)=1.e3
1939	qsurf_hist(ig,igcm_h2o_vap)=0.0
1940	endif
1941	endif
1942
1943	enddo
1944	endif
1945
1946	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1947	call physdem1("restartfi.nc",long,lati,nsoilmx,nq, &
1948	! ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf, &
1949	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1950	area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1951	cloudfrac,totcloudfrac,hice)
1952	endif
1953
1954	! -----------------------------------------------------------------
1955	! Saving statistics :
1956	! -----------------------------------------------------------------
1957	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1958	! which can later be used to make the statistic files of the run:
1959	! "stats") only possible in 3D runs !
1960
1961
1962	if (callstats) then
1963
1964	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1965	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1966	call wstats(ngrid,"fluxsurf_lw", &
1967	"Thermal IR radiative flux to surface","W.m-2",2, &
1968	fluxsurf_lw)
1969	! call wstats(ngrid,"fluxsurf_sw", &
1970	! "Solar radiative flux to surface","W.m-2",2, &
1971	! fluxsurf_sw_tot)
1972	call wstats(ngrid,"fluxtop_lw", &
1973	"Thermal IR radiative flux to space","W.m-2",2, &
1974	fluxtop_lw)
1975	! call wstats(ngrid,"fluxtop_sw", &
1976	! "Solar radiative flux to space","W.m-2",2, &
1977	! fluxtop_sw_tot)
1978	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1979	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1980	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1981	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1982	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1983
1984	if (tracer) then
1985	if (water) then
1986	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
1987	call wstats(ngrid,"vmr_h2ovapor", &
1988	"H2O vapour volume mixing ratio","mol/mol", &
1989	3,vmr)
1990	endif ! of if (water)
1991
1992	endif !tracer
1993
1994	if(lastcall) then
1995	write (,) "Writing stats..."
1996	call mkstats(ierr)
1997	endif
1998	endif !if callstats
1999
2000
2001	! ----------------------------------------------------------------------
2002	! output in netcdf file "DIAGFI", containing any variable for diagnostic
2003	! (output with period "ecritphy", set in "run.def")
2004	! ----------------------------------------------------------------------
2005	! writediagfi can also be called from any other subroutine for any variable.
2006	! but its preferable to keep all the calls in one place...
2007
2008	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2009	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2010
2011	if(1.eq.2)then
2012	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2013	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2014	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2015	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2016	call writediagfi(ngrid,'p','Pressure','Pa',3,pplay)
2017	endif
2018
2019	! Total energy balance diagnostics
2020	if(callrad.and.(.not.newtonian))then
2021	call writediagfi(ngrid,'ALB','Surface albedo',' ',2,albedo)
2022	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2023	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2024	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2025	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2026	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2027	!call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
2028	!call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
2029	!call writediagfi(ngrid,"vdifcdE","heat from vdifc surface","W m-2",2,vdifcdE)
2030	endif
2031
2032	! Temporary inclusions for heating diagnostics
2033	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2034	! call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2035	! call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2036	! call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2037
2038	! debugging
2039	!call writediagfi(ngrid,"vdifNC","H2O loss vdifc","kg m-2 s-1",2,vdifcncons)
2040	!call writediagfi(ngrid,"cadjNC","H2O loss convadj","kg m-2 s-1",2,cadjncons)
2041	call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2042	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2043	!print*,vdifcncons
2044	!call abort
2045
2046	! Output aerosols
2047	! call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
2048	! call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
2049	! call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
2050	! call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
2051
2052	! Output tracers
2053	if (tracer) then
2054	do iq=1,nq
2055	! call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2056	! call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2057	! 'kg m^-2',2,qsurf(1,iq) )
2058	call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2059	'kg m^-2',2,qsurf_hist(1,iq) )
2060	call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2061	'kg m^-2',2,qcol(1,iq) )
2062
2063	if(water)then
2064	call writediagfi(ngridmx,"H2Omaxcol","max. poss. H2O column","kg m^-2",2,H2Omaxcol)
2065	endif
2066
2067	if(watercond)then
2068	!call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2069	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2070	endif
2071
2072	if(waterrain)then
2073	call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
2074	call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
2075	endif
2076
2077	if(hydrology)then
2078	call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
2079	endif
2080
2081	if(ice_update)then
2082	call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
2083	call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
2084	endif
2085
2086	do ig=1,ngrid
2087	if(tau_col(ig).gt.1.e3)then
2088	print*,'WARNING: tau_col=',tau_col(ig)
2089	print*,'at ig=',ig,'in PHYSIQ'
2090	!call abort
2091	endif
2092	end do
2093
2094	call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2095
2096	enddo
2097	endif
2098
2099	! ----------------------------------------------------------
2100	! Output in netcdf file "diagsoil.nc" for subterranean
2101	! variables (output every "ecritphy", as for writediagfi)
2102	! ----------------------------------------------------------
2103
2104	! Write soil temperature
2105	!call writediagsoil(ngrid,"soiltemp","Soil temperature","K",3,tsoil)
2106
2107	else ! if(ngrid.eq.1)
2108
2109	! ----------------------------------------------------------------------
2110	! Output in grads file "g1d" (ONLY when using testphys1d)
2111	! ----------------------------------------------------------------------
2112
2113	if(countG1D.eq.saveG1D)then
2114
2115	call writeg1d(ngrid,1,tsurf,'tsurf','K')
2116	call writeg1d(ngrid,1,ps,'ps','Pa')
2117	call writeg1d(ngrid,nlayer,zt,'T','K')
2118	call writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg')
2119	call writeg1d(ngrid,nlayer,aerosol,'aerosol','SI')
2120	call writeg1d(ngrid,nlayer,reffrad,'reffrad','SI')
2121	call writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI')
2122	call writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI')
2123	call writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI')
2124
2125	! radiation balance
2126	call writeg1d(ngrid,1,ISR/totarea,'ISR','W m-2')
2127	call writeg1d(ngrid,1,ASR/totarea,'ASR','W m-2')
2128	call writeg1d(ngrid,1,OLR/totarea,'OLR','W m-2')
2129
2130	if(tracer) then
2131	do iq=1,nq
2132	call writeg1d(ngrid,1,qsurf(1,iq),trim(noms(iq))//'_s','kg m^-2')
2133	call writeg1d(ngrid,1,qcol(1,iq),trim(noms(iq))//'_c','kg m^-2')
2134	call writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
2135	end do
2136
2137	if(waterrain)then
2138	CALL writeg1d(ngrid,1,zdqsrain,'rainfall','kg m-2 s-1')
2139	CALL writeg1d(ngrid,1,zdqssnow,'snowfall','kg m-2 s-1')
2140	endif
2141
2142	end if
2143
2144	countG1D=1
2145	else
2146	countG1D=countG1D+1
2147	endif ! if time to save
2148
2149	endif ! if(ngrid.ne.1)
2150
2151	icount=icount+1
2152
2153	return
2154	end subroutine physiq

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

Download in other formats:

Original Format