Context Navigation

physiq.F90 @ 961

Last change on this file since 961 was 961, checked in by jleconte, 12 years ago

15/05/2013 == JL

correction in radiative scheme to enforce double precision
corrects calculation of ISR

Property svn:executable set to *

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

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

Download in other formats:

Original Format