Context Navigation

physiq.F90 @ 538

Last change on this file since 538 was 538, checked in by rwordsworth, 13 years ago
calc_cpp_mugaz changed to a check only in 3D
Property svn:executable set to ``*
File size: 76.7 KB

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

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

Download in other formats:

Original Format