Context Navigation

physiq.F90 @ 862

Last change on this file since 862 was 862, checked in by aslmd, 13 years ago

LMDZ.UNIVERSAL
LMDZ.GENERIC

Added calls to Frederic Hourdin's subroutines to output physical fields in parallel. This is under precompiling flags CPP_PARA.
This allows for LMDZ.UNIVERSAL users to use writediagfi with parallel computations.
These lines are not compiled by casual users of LMDZ.GENERIC (or users of LMDZ.UNIVERSAL in sequential mode).

The precompiling flag NOWRITEDIAGFI is obsolete and has been deleted.

NOTES:

A better cleaner method might be proposed later
Added subroutines

A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_ecrit.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iophys.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd.h
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_ini.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_fin.F90
The iotd_* subroutines are actually supposed to be put in bibio in LMDZ.COMMON
This will be done later once agreed in the team

Property svn:executable set to *

File size: 75.8 KB

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

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

Download in other formats:

Original Format