Context Navigation

physiq.F90 @ 858

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

Generic GCM:

Fixed sedimentation issue: ensure in callsedim that the correct radii are provided to newsedim and also that the updated temperature and tracer mixing ratios are used to compute sedimentation.
Updated the way aerosol radii are considered and used; routines in radii_mod (h2o_reffrad, co2_reffrad, etc.) only handle a single aerosol. The idea here is that these can be called from anywhere and that the caller doesn't need to have the full (naerkind size) array of aerosol radii.
cleanup (addition of intent(..) to routine arguments) in various routines

Property svn:executable set to *

File size: 75.6 KB

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

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

Download in other formats:

Original Format