Context Navigation

physiq.F90 @ 711

Last change on this file since 711 was 651, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

Property svn:executable set to *

File size: 69.5 KB

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

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

Download in other formats:

Original Format