Context Navigation

physiq.F90 @ 937

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

28/03/2013 == JL

optimization of optci and optcv routines. 15to 25% gain on these routines.

around 10% on the whole code with 1 scatterer.

No changes on output (byte to byte)

corrected bug in gray case in callcorrk.
added profiling option in makegcm_ifort. See the file for details

Property svn:executable set to *

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

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

Download in other formats:

Original Format