Context Navigation

physiq.F90 @ 848

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

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

Download in other formats:

Original Format