Context Navigation

physiq.F90 @ 486

Last change on this file since 486 was 486, checked in by rwordsworth, 13 years ago
Variable ice timestep added for ice evolution algorithm. Sedimentation improved to allow consistent dust (L. Kerber). Bug linked to allocatable matrices removed from rcm1d.F. Treatment of initial aerosol radii in callcorrk.F90 improved.
Property svn:executable set to ``*
File size: 77.5 KB

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

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

Download in other formats:

Original Format