Context Navigation

physiq.F90 @ 646

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

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

Download in other formats:

Original Format