Context Navigation

physiq.F90 @ 537

Last change on this file since 537 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

Property svn:executable set to *

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

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

Download in other formats:

Original Format