Context Navigation

physiq.F90 @ 601

Last change on this file since 601 was 600, checked in by jleconte, 13 years ago

Corrects the computation of planck function at the surface in sfluxi

so that its integral is equal to sigma Tsurf^4.

This ensure that no flux is lost due to:

-truncation of the planck function at high/low wavenumber
-numerical error during first spectral computation of the planck function
-discrepancy between Tsurf and NTS/NTfac in sfluxi

OLR now equal to LW net heating/cooling at equilibrium!
As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the

corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.

Property svn:executable set to *

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

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

Download in other formats:

Original Format