Context Navigation

physiq.F90 @ 314

Last change on this file since 314 was 305, checked in by rwordsworth, 14 years ago
Several new files added as part of the climate evolution model (main program kcm.F90). Some general cleanup in physiq.F90 and callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh scattering corrected.
Property svn:executable set to ``*
File size: 77.0 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
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
455	if (pday.ne.day_ini) then
456	write(,) "ERROR in physiq.F90:"
457	write(,) "bad synchronization between physics and dynamics"
458	write(,) "dynamics day: ",pday
459	write(,) "physics day: ",day_ini
460	stop
461	endif
462
463	write (,) 'In physiq day_ini =', day_ini
464
465	! Initialize albedo and orbital calculation
466	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
467	call surfini(ngrid,qsurf,albedo0)
468	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
469
470	do ig=1,ngrid
471	albedo(ig)=albedo0(ig)
472	enddo
473
474	if(tlocked)then
475	print*,'Planet is tidally locked at resonance n=',nres
476	print*,'Make sure you have the right rotation rate!!!'
477	endif
478
479	! initialize soil
480	! ~~~~~~~~~~~~~~~
481	if (callsoil) then
482	call soil(ngrid,nsoilmx,firstcall,inertiedat, &
483	ptimestep,tsurf,tsoil,capcal,fluxgrd)
484	else
485	print*,'WARNING! Thermal conduction in the soil turned off'
486	do ig=1,ngrid
487	capcal(ig)=1.e6
488	fluxgrd(ig)=0.
489	if(noradsurf)then
490	fluxgrd(ig)=10.0
491	endif
492	enddo
493	print*,'Flux from ground = ',fluxgrd,' W m^-2'
494	endif
495	icount=1
496
497	! decide whether to update ice at end of run
498	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
499	ice_update=.false.
500	if(sourceevol)then
501	open(128,file='num_run',form='formatted')
502	read(128,*) num_run
503	close(128)
504
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	! print*,'OCEAN EVERYWHERE FOR GJ581d RUNS'
515	do ig=1,ngridmx
516	rnat(ig)=1
517
518	! if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
519	if(inertiedat(ig,1).gt.1.E4)then
520	rnat(ig)=0
521	endif
522	enddo
523
524	print*,'WARNING! Surface type currently decided by surface inertia'
525	print*,'This should be improved e.g. in newstart.F'
526
527
528	! initialise surface history variable
529	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
530	do ig=1,ngridmx
531	do iq=1,nqmx
532	qsurf_hist(ig,iq)=qsurf(ig,iq)
533	enddo
534	enddo
535
536	! initialise variable for dynamical heating diagnostic
537	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
538	ztprevious(:,:)=pt(:,:)
539
540	! Set temperature just above condensation temperature (for Early Mars)
541	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
542	if(nearco2cond) then
543	write(,)' WARNING! Starting at Tcond+1K'
544	do l=1, nlayer
545	do ig=1,ngrid
546	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
547	-pt(ig,l)) / ptimestep
548	enddo
549	enddo
550	endif
551
552	if(meanOLR)then
553	! to record global radiative balance
554	call system('rm -f rad_bal.out')
555	! to record global mean/max/min temperatures
556	call system('rm -f tem_bal.out')
557	! to record global hydrological balance
558	call system('rm -f h2o_bal.out')
559	endif
560
561	watertest=.false.
562	if(water)then
563	! initialise variables for water cycle
564
565	qzero1D=0.0
566
567	watertest=.true. ! standard for now
568	if(ngrid.eq.1)then
569	qzero1D = 0.0
570	qsurf(1,igcm_h2o_vap) = qzero1D
571	do l=1, nlayer
572	pq(1,l,igcm_h2o_vap)=0.0
573	pq(1,l,igcm_h2o_ice)=0.0
574	enddo
575	endif
576
577	do ig=1,ngrid
578	qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
579	if(ice_update)then
580	ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
581	endif
582	enddo
583
584	endif
585	call su_watercycle ! even if we don't have a water cycle, we might
586	! need epsi for the wvp definitions in callcorrk.F
587
588	endif ! (end of "if firstcall")
589
590	! ---------------------------------------------------
591	! 1.2 Initializations done at every physical timestep:
592	! ---------------------------------------------------
593
594	if (ngrid.NE.ngridmx) then
595	print*,'STOP in PHYSIQ'
596	print*,'Probleme de dimensions :'
597	print*,'ngrid = ',ngrid
598	print*,'ngridmx = ',ngridmx
599	stop
600	endif
601
602	! Initialize various variables
603	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
604
605	pdu(:,:) = 0.0
606	pdv(:,:) = 0.0
607	! if ( (.not.nearco2cond).and.(.not.firstcall) ) then
608	if ( .not.nearco2cond ) then
609	pdt(:,:) = 0.0
610	endif ! this was the source of an evil bug...
611	pdq(:,:,:) = 0.0
612	pdpsrf(:) = 0.0
613	zflubid(:) = 0.0
614	zdtsurf(:) = 0.0
615	dqsurf(:,:) = 0.0
616
617	zday=pday+ptime ! compute time, in sols (and fraction thereof)
618
619	! Compute Stellar Longitude (Ls)
620	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
621	if (season) then
622	call stellarlong(zday,zls)
623	else
624	call stellarlong(float(day_ini),zls)
625	end if
626
627	! Compute geopotential between layers
628	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
629
630	do l=1,nlayer
631	do ig=1,ngrid
632	zzlay(ig,l)=pphi(ig,l)/g
633	enddo
634	enddo
635	do ig=1,ngrid
636	zzlev(ig,1)=0.
637	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
638	enddo
639	do l=2,nlayer
640	do ig=1,ngrid
641	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
642	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
643	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
644	enddo
645	enddo
646	! Potential temperature calculation may not be the same in physiq and dynamic...
647
648	! Compute potential temperature
649	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
650
651	do l=1,nlayer
652	do ig=1,ngrid
653	call calc_cpp3d(cppNI,rcpNI,pt(ig,l),pplay(ig,l))
654	cpp3D(ig,l) = cppNI
655	rcp3D(ig,l) = rcpNI
656	enddo
657	enddo
658
659	do l=1,nlayer
660	do ig=1,ngrid
661	zh(ig,l) = pt(ig,l)
662	zpopskNI(ig,l) = pt(ig,l)/zh(ig,l)
663	zpopsk(ig,l) = pt(ig,l)/zh(ig,l)
664	! we're only after zpopskNI here, not zh
665	! zh is calculated seperately before both vdifc and convadj
666	enddo
667	enddo
668
669	do l=1,nlayer
670	do ig=1,ngrid
671	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
672	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
673	enddo
674	enddo
675
676	!-----------------------------------------------------------------------
677	! 2. Compute radiative tendencies
678	!-----------------------------------------------------------------------
679
680	if (callrad) then
681	if( mod(icount-1,iradia).eq.0) then
682
683	! Compute local stellar zenith angles
684	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
685	call orbite(zls,dist_star,declin)
686
687	if (tlocked) then
688	ztim1=SIN(declin)
689	ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
690	ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
691
692	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
693	ztim1,ztim2,ztim3,mu0,fract)
694
695	elseif (diurnal) then
696	ztim1=SIN(declin)
697	ztim2=COS(declin)COS(2.pi*(zday-.5))
698	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
699
700	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
701	ztim1,ztim2,ztim3,mu0,fract)
702
703	else
704
705	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
706	! WARNING: this function appears not to work in 1D
707
708	endif
709
710	if (corrk) then
711
712	! a) Call correlated-k radiative transfer scheme
713	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
714
715	if(kastprof)then
716	print*,'kastprof should not = true here'
717	call abort
718	endif
719	muvar(:,:)=0.0 ! only used for climate evolution for now
720
721	! Option to call scheme once for clear regions, once for cloudy regions (BC)
722	if(CLFvarying)then
723
724	clearsky=.true.
725	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
726	albedo,emis,mu0,pplev,pplay,pt, &
727	tsurf,fract,dist_star,aerosol,cpp3D,muvar, &
728	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
729	fluxabs_sw1,fluxtop_dn,reffrad,tau_col1, &
730	cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
731
732	clearsky=.false.
733	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
734	albedo,emis,mu0,pplev,pplay,pt, &
735	tsurf,fract,dist_star,aerosol,cpp3D,muvar, &
736	zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, &
737	fluxabs_sw2,fluxtop_dn,reffrad,tau_col2, &
738	cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
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_lw2(ig)
746	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw2(ig)
747	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw2(ig)
748	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw2(ig)
749	tau_col(ig) = ntftau_col1(ig) + tftau_col2(ig)
750
751	do l=1,nlayer
752	zdtlw(ig,l) = ntfzdtlw1(ig,l) + tfzdtlw2(ig,l)
753	zdtsw(ig,l) = ntfzdtsw1(ig,l) + tfzdtsw2(ig,l)
754	enddo
755	enddo
756
757	! Otherwise standard callcorrk
758	else
759
760	clearsky=.false.
761	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
762	albedo,emis,mu0,pplev,pplay,pt, &
763	tsurf,fract,dist_star,aerosol,cpp3D,muvar, &
764	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
765	fluxabs_sw,fluxtop_dn,reffrad,tau_col, &
766	cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
767
768	endif
769
770	! Radiative flux from the sky absorbed by the surface (W.m-2)
771	GSR=0.0
772	do ig=1,ngrid
773	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig) &
774	+fluxsurf_sw(ig)*(1.-albedo(ig))
775
776	if(noradsurf)then ! no lower surface; SW flux just disappears
777	GSR = GSR + fluxsurf_sw(ig)*area(ig)
778	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
779	endif
780
781	enddo
782	if(noradsurf)then
783	print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
784	endif
785
786	! Net atmospheric radiative heating rate (K.s-1)
787	do l=1,nlayer
788	do ig=1,ngrid
789	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
790	enddo
791	enddo
792
793
794	elseif(newtonian)then
795
796	! b) Call Newtonian cooling scheme
797	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
798	call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
799
800	do ig=1,ngrid
801	zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
802	enddo
803	! e.g. surface becomes proxy for 1st atmospheric layer ?
804
805	else
806
807	! c) Atmosphere has no radiative effect
808	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
809	do ig=1,ngrid
810	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_star**2
811	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
812	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
813	endif
814	fluxsurf_sw(ig) = fluxtop_dn(ig)
815	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
816	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
817	enddo ! radiation skips the atmosphere entirely
818
819	do l=1,nlayer
820	do ig=1,ngrid
821	dtrad(ig,l)=0.0
822	enddo
823	enddo ! hence no atmospheric radiative heating
824
825	endif ! if corrk
826
827	endif ! of if(mod(icount-1,iradia).eq.0)
828
829
830	! Transformation of the radiative tendencies
831	! ------------------------------------------
832
833	do ig=1,ngrid
834	zplanck(ig)=tsurf(ig)*tsurf(ig)
835	zplanck(ig)=emis(ig)stephanzplanck(ig)*zplanck(ig)
836	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
837	enddo
838
839	do l=1,nlayer
840	do ig=1,ngrid
841	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
842	enddo
843	enddo
844
845	!-------------------------
846	! test energy conservation
847	if(enertest)then
848	dEtotSW = 0.0
849	dEtotLW = 0.0
850	dEtotsSW = 0.0
851	dEtotsLW = 0.0
852	do ig = 1, ngrid
853	do l = 1, nlayer
854	masse = (pplev(ig,l) - pplev(ig,l+1))/g
855	dEtotSW = dEtotSW + cpp3D(ig,l)massezdtsw(ig,l)*area(ig)
856	dEtotLW = dEtotLW + cpp3D(ig,l)massezdtlw(ig,l)*area(ig)
857	enddo
858	dEtotsSW = dEtotsSW + fluxsurf_sw(ig)(1.-albedo(ig))area(ig)
859	dEtotsLW = dEtotsLW + emis(ig)fluxsurf_lw(ig)area(ig)-zplanck(ig)*area(ig)
860	enddo
861	dEtotSW = dEtotSW/totarea
862	dEtotLW = dEtotLW/totarea
863	dEtotsSW = dEtotsSW/totarea
864	dEtotsLW = dEtotsLW/totarea
865	print*,'---------------------------------------------------------------'
866	print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2'
867	print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2'
868	print*,'In corrk SW surface energy change =',dEtotsSW,' W m-2'
869	print*,'In corrk LW surface energy change =',dEtotsLW,' W m-2'
870	endif
871	!-------------------------
872
873	endif ! of if (callrad)
874
875	!-----------------------------------------------------------------------
876	! 4. Vertical diffusion (turbulent mixing):
877	! -----------------------------------------
878
879	if (calldifv) then
880
881	do ig=1,ngrid
882	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
883	enddo
884
885	zdum1(:,:)=0.0
886	zdum2(:,:)=0.0
887	do l=1,nlayer
888	do ig=1,ngrid
889	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
890	enddo
891	enddo
892
893	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
894	ptimestep,capcal,lwrite, &
895	pplay,pplev,zzlay,zzlev,z0, &
896	pu,pv,zh,pq,tsurf,emis,qsurf, &
897	zdum1,zdum2,zdh,pdq,zflubid, &
898	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
899	zdqdif,zdqsdif,lastcall)
900
901	do l=1,nlayer
902	do ig=1,ngrid
903	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
904	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
905	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
906	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
907	enddo
908	enddo
909
910	do ig=1,ngrid
911	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
912	enddo
913
914	if (tracer) then
915	do iq=1, nq
916	do l=1,nlayer
917	do ig=1,ngrid
918	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
919	enddo
920	enddo
921	enddo
922	do iq=1, nq
923	do ig=1,ngrid
924	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
925	enddo
926	enddo
927
928	end if ! of if (tracer)
929
930	!-------------------------
931	! test energy conservation
932	if(enertest)then
933	dEtot=0.0
934	dEtots=0.0
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	dEtots = dEtots + capcal(ig)zdtsdif(ig)area(ig)
949	vdifcdE(ig) = vdifcdE(ig) + capcal(ig)*zdtsdif(ig)
950	enddo
951	dEtot = dEtot/totarea
952	dEtots = dEtots/totarea
953	print*,'In difv atmospheric energy change =',dEtot,' W m-2'
954	print*,'In difv surface energy change =',dEtots,' W m-2'
955	!print*,'Note we ignore the wind change...'
956	! and latent heat for that matter...
957	endif
958	!-------------------------
959
960	!-------------------------
961	! test water conservation
962	if(watertest.and.water)then
963	dWtot=0.0
964	dWtots=0.0
965	nconsMAX=0.0
966	do ig = 1, ngrid
967
968	vdifcncons(ig)=0.0
969	do l = 1, nlayer
970	masse = (pplev(ig,l) - pplev(ig,l+1))/g
971
972	iq = igcm_h2o_vap
973	dWtot = dWtot + massezdqdif(ig,l,iq)ptimestep*area(ig)
974	vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
975
976	iq = igcm_h2o_ice
977	dWtot = dWtot + massezdqdif(ig,l,iq)ptimestep*area(ig)
978	vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
979
980	enddo
981
982	iq = igcm_h2o_vap
983	dWtots = dWtots + zdqsdif(ig,iq)ptimesteparea(ig)
984	vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
985
986	iq = igcm_h2o_ice
987	dWtots = dWtots + zdqsdif(ig,iq)ptimesteparea(ig)
988	vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
989
990	if(vdifcncons(ig).gt.nconsMAX)then
991	nconsMAX=vdifcncons(ig)
992	endif
993
994	enddo
995
996	dWtot = dWtot/totarea
997	dWtots = dWtots/totarea
998	print*,'---------------------------------------------------------------'
999	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
1000	print*,'In difv surface water change =',dWtots,' kg m-2'
1001	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
1002	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
1003
1004
1005	endif
1006	!-------------------------
1007
1008	else
1009
1010	if(.not.newtonian)then
1011
1012	do ig=1,ngrid
1013	zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
1014	enddo
1015
1016	endif
1017
1018	endif ! of if (calldifv)
1019
1020
1021	!-----------------------------------------------------------------------
1022	! 5. Dry convective adjustment:
1023	! -----------------------------
1024
1025	if(calladj) then
1026
1027	do l=1,nlayer
1028	do ig=1,ngrid
1029	if(nonideal)then
1030	zdh(ig,l) = pdt(ig,l)/zpopskNI(ig,l)
1031	else
1032	zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
1033	endif
1034	enddo
1035	enddo
1036	zduadj(:,:)=0.0
1037	zdvadj(:,:)=0.0
1038	zdhadj(:,:)=0.0
1039
1040	if(nonideal)then
1041	print*,'Nonideal is not used at the moment'
1042	call abort
1043	call convadj(ngrid,nlayer,nq,ptimestep, &
1044	pplay,pplev,zpopskNI, &
1045	pu,pv,zh,pq, &
1046	pdu,pdv,zdh,pdq, &
1047	zduadj,zdvadj,zdhadj, &
1048	zdqadj)
1049	else
1050
1051	call convadj(ngrid,nlayer,nq,ptimestep, &
1052	pplay,pplev,zpopsk, &
1053	pu,pv,zh,pq, &
1054	pdu,pdv,zdh,pdq, &
1055	zduadj,zdvadj,zdhadj, &
1056	zdqadj)
1057
1058	endif
1059
1060	do l=1,nlayer
1061	do ig=1,ngrid
1062	pdu(ig,l) = pdu(ig,l) + zduadj(ig,l)
1063	pdv(ig,l) = pdv(ig,l) + zdvadj(ig,l)
1064	if(nonideal)then
1065	pdt(ig,l) = pdt(ig,l) + zdhadj(ig,l)*zpopskNI(ig,l)
1066	zdtadj(ig,l) = zdhadj(ig,l)*zpopskNI(ig,l) ! for diagnostic only
1067	else
1068	pdt(ig,l) = pdt(ig,l) + zdhadj(ig,l)*zpopsk(ig,l)
1069	zdtadj(ig,l) = zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1070	endif
1071	enddo
1072	enddo
1073
1074	if(tracer) then
1075	do iq=1, nq
1076	do l=1,nlayer
1077	do ig=1,ngrid
1078	pdq(ig,l,iq) = pdq(ig,l,iq) + zdqadj(ig,l,iq)
1079	enddo
1080	enddo
1081	enddo
1082	end if
1083
1084	!-------------------------
1085	! test energy conservation
1086	if(enertest)then
1087	dEtot=0.0
1088	do ig = 1, ngrid
1089	do l = 1, nlayer
1090	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1091	dEtot = dEtot + cpp3D(ig,l)massezdtadj(ig,l)*area(ig)
1092	enddo
1093	enddo
1094	dEtot=dEtot/totarea
1095	print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1096	endif
1097	!-------------------------
1098
1099	!-------------------------
1100	! test water conservation
1101	if(watertest)then
1102	dWtot=0.0
1103	do ig = 1, ngrid
1104	cadjncons(ig)=0.0
1105	do l = 1, nlayer
1106	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1107
1108	iq = igcm_h2o_vap
1109	dWtot = dWtot + massezdqadj(ig,l,iq)ptimestep*area(ig)
1110	cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
1111
1112	iq = igcm_h2o_ice
1113	dWtot = dWtot + massezdqadj(ig,l,iq)ptimestep*area(ig)
1114	cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
1115
1116	enddo
1117
1118	enddo
1119	dWtot=dWtot/totarea
1120	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
1121	endif
1122	!-------------------------
1123
1124	endif ! of if(calladj)
1125
1126	!-----------------------------------------------------------------------
1127	! 6. Carbon dioxide condensation-sublimation:
1128	! -------------------------------------------
1129
1130	if (co2cond) then
1131	if (.not.tracer) then
1132	print*,'We need a CO2 ice tracer to condense CO2'
1133	call abort
1134	endif
1135
1136	call condense_cloud(ngrid,nlayer,nq,ptimestep, &
1137	capcal,pplay,pplev,tsurf,pt, &
1138	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1139	qsurf(1,igcm_co2_ice),albedo,emis, &
1140	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1141	zdqc,reffrad,cpp3D)
1142
1143	do l=1,nlayer
1144	do ig=1,ngrid
1145	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1146	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1147	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1148	enddo
1149	enddo
1150	do ig=1,ngrid
1151	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1152	enddo
1153
1154	do iq=1,nq ! should use new notation here !
1155	do l=1,nlayer
1156	do ig=1,ngrid
1157	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1158	enddo
1159	enddo
1160	enddo
1161	! Note: we do not add surface co2ice tendency
1162	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
1163
1164	!-------------------------
1165	! test energy conservation
1166	if(enertest)then
1167	dEtot=0.0
1168	dEtots=0.0
1169	do ig = 1, ngrid
1170	do l = 1, nlayer
1171	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1172	dEtot = dEtot + cpp3D(ig,l)massezdtc(ig,l)*area(ig)
1173	enddo
1174	dEtots = dEtots + capcal(ig)zdtsurfc(ig)area(ig)
1175	enddo
1176	dEtot=dEtot/totarea
1177	dEtots=dEtots/totarea
1178	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
1179	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
1180	endif
1181	!-------------------------
1182
1183	endif ! of if (co2cond)
1184
1185
1186	!-----------------------------------------------------------------------
1187	! 7. Specific parameterizations for tracers
1188	! -----------------------------------------
1189
1190	if (tracer) then
1191
1192	! 7a. Water and ice
1193	! ---------------
1194	if (water) then
1195
1196	! ----------------------------------------
1197	! Water ice condensation in the atmosphere
1198	! ----------------------------------------
1199	if(watercond)then
1200
1201	if(RLVTT.gt.1.e-8)then
1202
1203	! Re-evaporate cloud water/ice
1204	call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1205	DO l = 1, nlayer
1206	DO ig = 1, ngrid
1207	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1208	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1209	pdt(ig,l) = pdt(ig,l)+dtevap(ig,l)
1210	enddo
1211	enddo
1212
1213	call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1214	do l=1,nlayer
1215	do ig=1,ngrid
1216	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
1217	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
1218	pdt(ig,l) = pdt(ig,l)+dtmoist(ig,l)
1219	enddo
1220	enddo
1221
1222	!-------------------------
1223	! test energy conservation
1224	if(enertest)then
1225	dEtot=0.0
1226	madjdE(:)=0.0
1227	do ig = 1, ngrid
1228	do l = 1, nlayer
1229	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1230	dEtot = dEtot + cpp3D(ig,l)masse(dtmoist(ig,l)+dtevap(ig,l))*area(ig)
1231	madjdE(ig) = madjdE(ig) + cpp3D(ig,l)masse(dtmoist(ig,l)+dtevap(ig,l))
1232	enddo
1233	enddo
1234	dEtot=dEtot/totarea
1235	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
1236	endif
1237	!-------------------------
1238
1239	!-------------------------
1240	! test water conservation
1241	if(watertest)then
1242	dWtot=0.0
1243	do ig = 1, ngrid
1244	do iq = 1 , nq
1245	do l = 1, nlayer
1246	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1247	dWtot = dWtot + massedqmoist(ig,l,igcm_h2o_vap)area(ig)*ptimestep
1248	dWtot = dWtot + massedqmoist(ig,l,igcm_h2o_ice)area(ig)*ptimestep
1249	enddo
1250	enddo
1251	enddo
1252	dWtot=dWtot/totarea
1253	print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
1254	endif
1255	!-------------------------
1256
1257
1258	endif
1259
1260
1261	! Re-evaporate cloud water/ice
1262	call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1263	do l = 1, nlayer
1264	do ig = 1, ngrid
1265	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1266	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1267	pdt(ig,l) = pdt(ig,l)+dtevap(ig,l)
1268	enddo
1269	enddo ! note: we use qevap but not tevap in largescale/moistadj
1270	! otherwise is a big mess
1271
1272	call largescale(ptimestep,pplev,pplay,pt,qevap, & ! a bug was here!
1273	pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
1274	do l=1,nlayer
1275	do ig=1,ngrid
1276	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
1277	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
1278	pdt(ig,l) = pdt(ig,l)+dtlscale(ig,l)
1279
1280	if(.not.aerofixed)then
1281	reffrad(ig,l,2)=reffH2O(ig,l)
1282	endif
1283
1284	enddo
1285	enddo
1286
1287	!-------------------------
1288	! test energy conservation
1289	if(enertest)then
1290	dEtot=0.0
1291	lscaledE(:)=0.0
1292	do ig = 1, ngrid
1293	do l = 1, nlayer
1294	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1295	dEtot = dEtot + cpp3D(ig,l)masse(dtlscale(ig,l)+dtevap(ig,l))*area(ig)
1296	lscaledE(ig) = lscaledE(ig) + cpp3D(ig,l)masse(dtlscale(ig,l)+dtevap(ig,l))
1297	enddo
1298	enddo
1299	dEtot=dEtot/totarea
1300	print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
1301	endif
1302	!-------------------------
1303
1304	!-------------------------
1305	! test water conservation
1306	if(watertest)then
1307	dWtot=0.0
1308	do ig = 1, ngrid
1309	do iq = 1 , nq
1310	do l = 1, nlayer
1311	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1312	dWtot = dWtot + massedqvaplscale(ig,l)area(ig)*ptimestep
1313	dWtot = dWtot + massedqcldlscale(ig,l)area(ig)*ptimestep
1314	enddo
1315	enddo
1316	enddo
1317	dWtot=dWtot/totarea
1318	print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
1319	endif
1320	!-------------------------
1321
1322	! compute cloud fraction
1323	do l = 1, nlayer
1324	do ig = 1,ngrid
1325	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1326	enddo
1327	enddo
1328
1329	! compute total cloud fraction in column
1330	call totalcloudfrac(cloudfrac,totcloudfrac)
1331
1332	endif ! of if (watercondense)
1333
1334
1335	! --------------------------------
1336	! Water ice / liquid precipitation
1337	! --------------------------------
1338	if(waterrain)then
1339
1340	zdqrain(:,:,:) = 0.0
1341	zdqsrain(:) = 0.0
1342	zdqssnow(:) = 0.0
1343
1344	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
1345	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
1346
1347	do l=1,nlayer
1348	do ig=1,ngrid
1349	pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
1350	pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
1351	pdt(ig,l) = pdt(ig,l)+zdtrain(ig,l)
1352	enddo
1353	enddo
1354
1355	do ig=1,ngrid
1356	dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
1357	dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
1358	rainout(ig) = zdqsrain(ig)+zdqssnow(ig) ! diagnostic
1359	enddo
1360
1361	!-------------------------
1362	! test energy conservation
1363	if(enertest)then
1364	dEtot=0.0
1365	do ig = 1, ngrid
1366	do l = 1, nlayer
1367	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1368	dEtot = dEtot + cpp3D(ig,l)massezdtrain(ig,l)*area(ig)
1369	enddo
1370	enddo
1371	dEtot=dEtot/totarea
1372	print*,'In rain atmospheric energy change =',dEtot,' W m-2'
1373	endif
1374	!-------------------------
1375
1376
1377	!-------------------------
1378	! test energy conservation
1379	if(enertest)then
1380	dEtot=0.0
1381	do ig = 1, ngrid
1382	do l = 1, nlayer
1383	masse = (pplev(ig,l) - pplev(ig,l+1))/g
1384	dEtot = dEtot + cpp3D(ig,l)massezdtrain(ig,l)*area(ig)
1385	enddo
1386	enddo
1387	dEtot=dEtot/totarea
1388	print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
1389
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	dItot = dItot + massezdqrain(ig,l,igcm_h2o_ice)area(ig)
1395	dVtot = dVtot + massezdqrain(ig,l,igcm_h2o_vap)area(ig)
1396	enddo
1397	dItot = dItot + zdqssnow(ig)*area(ig)
1398	dVtot = dVtot + zdqsrain(ig)*area(ig)
1399	enddo
1400	dEtot=(dItotRLVTT/cpp + dVtotRLVTT/cpp)/totarea
1401	print,'In rain dItot =',dItotRLVTT/(cpp*totarea),' W m-2'
1402	print,'In rain dVtot =',dVtotRLVTT/(cpp*totarea),' W m-2'
1403	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
1404	endif
1405	!-------------------------
1406
1407	!-------------------------
1408	! test water conservation
1409	if(watertest)then
1410	dWtot=0.0
1411	dWtots=0.0
1412	do ig = 1, ngrid
1413	!do iq = 1 , nq
1414	do l = 1, nlayer
1415	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1416	dWtot = dWtot + massezdqrain(ig,l,igcm_h2o_vap)area(ig)*ptimestep
1417	dWtot = dWtot + massezdqrain(ig,l,igcm_h2o_ice)area(ig)*ptimestep
1418	enddo
1419	!enddo
1420	dWtots = dWtots + (zdqsrain(ig)+zdqssnow(ig))area(ig)ptimestep
1421	enddo
1422	dWtot=dWtot/totarea
1423	dWtots=dWtots/totarea
1424	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
1425	print*,'In rain surface water change =',dWtots,' kg m-2'
1426	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
1427	endif
1428	!-------------------------
1429
1430	end if ! of if (waterrain)
1431	end if ! of if (water)
1432
1433
1434	! 7c. Aerosol particles
1435	! -------------------
1436	! -------------
1437	! Sedimentation
1438	! -------------
1439	if (sedimentation) then
1440	zdqsed(:,:,:) = 0.0
1441	zdqssed(:,:) = 0.0
1442
1443
1444	!-------------------------
1445	! find qtot
1446	if(watertest)then
1447	dWtot=0.0
1448	dWtots=0.0
1449	iq=3
1450	do ig = 1, ngrid
1451	do l = 1, nlayer
1452	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1453	dWtot = dWtot + massepq(ig,l,iq)area(ig)*ptimestep
1454	dWtots = dWtots + massepdq(ig,l,iq)area(ig)*ptimestep
1455	enddo
1456	enddo
1457	dWtot=dWtot/totarea
1458	dWtots=dWtots/totarea
1459	print*,'Before sedim pq =',dWtot,' kg m-2'
1460	print*,'Before sedim pdq =',dWtots,' kg m-2'
1461	endif
1462	!-------------------------
1463
1464	call callsedim(ngrid,nlayer,ptimestep, &
1465	pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
1466
1467	!-------------------------
1468	! find qtot
1469	if(watertest)then
1470	dWtot=0.0
1471	dWtots=0.0
1472	iq=3
1473	do ig = 1, ngrid
1474	do l = 1, nlayer
1475	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1476	dWtot = dWtot + massepq(ig,l,iq)area(ig)*ptimestep
1477	dWtots = dWtots + massepdq(ig,l,iq)area(ig)*ptimestep
1478	enddo
1479	enddo
1480	dWtot=dWtot/totarea
1481	dWtots=dWtots/totarea
1482	print*,'After sedim pq =',dWtot,' kg m-2'
1483	print*,'After sedim pdq =',dWtots,' kg m-2'
1484	endif
1485	!-------------------------
1486
1487	do iq=1,nq
1488	! for now, we only allow H2O ice to sediment
1489	! and as in rain.F90, whether it falls as rain or snow depends
1490	! only on the surface temperature
1491	do ig=1,ngrid
1492	do l=1,nlayer
1493	pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
1494	enddo
1495	dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
1496	enddo
1497	enddo
1498
1499	!-------------------------
1500	! test water conservation
1501	if(watertest)then
1502	dWtot=0.0
1503	dWtots=0.0
1504	do iq=1,nq
1505	do ig = 1, ngrid
1506	do l = 1, nlayer
1507	masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
1508	dWtot = dWtot + massezdqsed(ig,l,iq)area(ig)*ptimestep
1509	enddo
1510	dWtots = dWtots + zdqssed(ig,iq)area(ig)ptimestep
1511	enddo
1512	enddo
1513	dWtot=dWtot/totarea
1514	dWtots=dWtots/totarea
1515	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
1516	print*,'In sedim surface ice change =',dWtots,' kg m-2'
1517	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
1518	endif
1519	!-------------------------
1520
1521	end if ! of if (sedimentation)
1522
1523
1524	! 7d. Updates
1525	! ---------
1526
1527	! ---------------------------------
1528	! Updating tracer budget on surface
1529	! ---------------------------------
1530
1531	if(hydrology)then
1532
1533	call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
1534	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
1535	! note: for now, also changes albedo in the subroutine
1536
1537	do ig=1,ngrid
1538	zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
1539	do iq=1,nq
1540	qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
1541	enddo
1542	enddo
1543	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
1544
1545	!-------------------------
1546	! test energy conservation
1547	if(enertest)then
1548	dEtot=0.0
1549	do ig = 1, ngrid
1550	dEtots = dEtots + capcal(ig)zdtsurf_hyd(ig)area(ig)
1551	enddo
1552	dEtot=dEtot/totarea
1553	print*,'In hydrol atmospheric energy change =',dEtot,' W m-2'
1554	endif
1555	!-------------------------
1556
1557	!-------------------------
1558	! test water conservation
1559	if(watertest)then
1560	dWtots=0.0
1561	do ig = 1, ngrid
1562	dWtots = dWtots + dqs_hyd(ig,igcm_h2o_ice)area(ig)ptimestep
1563	enddo
1564	dWtots=dWtots/totarea
1565	print*,'In hydrol surface ice change =',dWtots,' kg m-2'
1566	dWtots=0.0
1567	do ig = 1, ngrid
1568	dWtots = dWtots + dqs_hyd(ig,igcm_h2o_vap)area(ig)ptimestep
1569	enddo
1570	dWtots=dWtots/totarea
1571	print*,'In hydrol surface water change =',dWtots,' kg m-2'
1572	print*,'---------------------------------------------------------------'
1573	endif
1574	!-------------------------
1575
1576	ELSE ! of if (hydrology)
1577
1578	do iq=1,nq
1579	do ig=1,ngrid
1580	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1581	enddo
1582	enddo
1583
1584	END IF ! of if (hydrology)
1585
1586	! Add qsurf to qsurf_hist, which is what we save in
1587	! diagfi.nc etc. At the same time, we set the water
1588	! content of ocean gridpoints back to zero, in order
1589	! to avoid rounding errors in vdifc, rain
1590	do ig = 1, ngrid
1591	do iq = 1, nq
1592	if(iq.eq.igcm_h2o_vap .and. rnat(ig).eq.0)then ! if liquid water and terrain = ocean
1593	qsurf_hist(ig,iq) = qsurf(ig,iq)
1594	!qsurf(ig,iq) = qcol(ig,iq)
1595	! the value of qsurf we choose here makes NO DIFFERENCE TO ANYTHING AT ALL
1596	else
1597	qsurf_hist(ig,iq) = qsurf(ig,iq)
1598	endif
1599	enddo
1600	enddo
1601
1602	if(ice_update)then
1603	do ig = 1, ngrid
1604	ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
1605	enddo
1606	endif
1607
1608	endif ! of if (tracer)
1609
1610
1611	! print*,'we got this far 1'
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	DYN = DYN + area(ig)*fluxdyn(ig)
1739
1740	if(fluxtop_dn(ig).lt.0.0)then
1741	print*,'fluxtop_dn has gone crazy'
1742	print*,'fluxtop_dn=',fluxtop_dn(ig)
1743	print*,'tau_col=',tau_col(ig)
1744	print*,'aerosol=',aerosol(ig,:,:)
1745	print*,'temp= ',pt(ig,:)
1746	print*,'pplay= ',pplay(ig,:)
1747	call abort
1748	endif
1749	end do
1750
1751	if(ngridmx.eq.1)then
1752	DYN=0.0
1753	endif
1754
1755	print*,' ISR ASR OLR GND DYN [W m^-2]'
1756	print*, ISR/totarea,ASR/totarea,OLR/totarea,GND/totarea,DYN/totarea
1757
1758	if(enertest)then
1759	print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
1760	endif
1761
1762	if(meanOLR)then
1763	if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
1764	! to record global radiative balance
1765	open(92,file="rad_bal.out",form='formatted',access='append')
1766	write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
1767	close(92)
1768	open(93,file="tem_bal.out",form='formatted',access='append')
1769	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1770	close(93)
1771	endif
1772	endif
1773
1774	endif
1775
1776	! ------------------------------------------------------------------
1777	! Diagnostic to test radiative-convective timescales in code
1778	! ------------------------------------------------------------------
1779	if(testradtimes)then
1780	open(38,file="tau_phys.out",form='formatted',access='append')
1781	ig=1
1782	do l=1,nlayer
1783	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1784	enddo
1785	close(38)
1786	print*,'As testradtimes enabled, exiting physics on first call'
1787	call abort
1788	endif
1789
1790	! ---------------------------------------------------------
1791	! Compute column amounts (kg m-2) if tracers are enabled
1792	! ---------------------------------------------------------
1793	if(tracer)then
1794	qcol(:,:)=0.0
1795	do iq=1,nq
1796	do ig=1,ngrid
1797	do l=1,nlayer
1798	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) * &
1799	(pplev(ig,l) - pplev(ig,l+1)) / g
1800	enddo
1801	enddo
1802	enddo
1803
1804	! not generalised for arbitrary aerosols yet!!!
1805	reffcol(:,:)=0.0
1806	do ig=1,ngrid
1807	do l=1,nlayer
1808	if(co2cond)then
1809	reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &
1810	reffrad(ig,l,1) * &
1811	(pplev(ig,l) - pplev(ig,l+1)) / g
1812	endif
1813	if(water)then
1814	reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &
1815	reffrad(ig,l,2) * &
1816	(pplev(ig,l) - pplev(ig,l+1)) / g
1817	endif
1818	enddo
1819	enddo
1820
1821	endif
1822
1823	! ---------------------------------------------------------
1824	! Test for water conservation if water is enabled
1825	! ---------------------------------------------------------
1826
1827	if(water)then
1828
1829	icesrf = 0.0
1830	liqsrf = 0.0
1831	icecol = 0.0
1832	vapcol = 0.0
1833
1834	h2otot = 0.0
1835	do ig=1,ngrid
1836
1837	icesrf = icesrf + area(ig)*qsurf_hist(ig,igcm_h2o_ice)
1838	liqsrf = liqsrf + area(ig)*qsurf_hist(ig,igcm_h2o_vap)
1839	icecol = icecol + area(ig)*qcol(ig,igcm_h2o_ice)
1840	vapcol = vapcol + area(ig)*qcol(ig,igcm_h2o_vap)
1841
1842	h2otot = h2otot + area(ig)* &
1843	(qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap) &
1844	+qsurf_hist(ig,igcm_h2o_ice)+qsurf_hist(ig,igcm_h2o_vap))
1845	end do
1846
1847	print*,' Total water amount [kg m^-2]: ',h2otot/totarea
1848	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
1849	print*, icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
1850
1851	if(meanOLR)then
1852	if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
1853	! to record global water balance
1854	open(98,file="h2o_bal.out",form='formatted',access='append')
1855	write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
1856	close(98)
1857	endif
1858	endif
1859
1860	endif
1861
1862	! ---------------------------------------------------------
1863	! Calculate RH for diagnostic if water is enabled
1864	! ---------------------------------------------------------
1865
1866	if(water)then
1867	do l = 1, nlayer
1868	do ig = 1, ngrid
1869	call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
1870	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
1871	enddo
1872	enddo
1873
1874	! compute maximum possible H2O column amount (100% saturation)
1875	do ig=1,ngrid
1876	H2Omaxcol(ig)=0.0
1877	do l=1,nlayer
1878	H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) * &
1879	(pplev(ig,l) - pplev(ig,l+1))/g
1880	enddo
1881	enddo
1882
1883	endif
1884
1885
1886	if (ngrid.ne.1) then
1887	print*,''
1888	print,'--> Ls =',zls180./pi
1889	! -------------------------------------------------------------------
1890	! Writing NetCDF file "RESTARTFI" at the end of the run
1891	! -------------------------------------------------------------------
1892	! Note: 'restartfi' is stored just before dynamics are stored
1893	! in 'restart'. Between now and the writting of 'restart',
1894	! there will have been the itau=itau+1 instruction and
1895	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1896	! thus we store for time=time+dtvr
1897
1898	if(lastcall) then
1899	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1900
1901
1902	! Update surface ice distribution to iterate to steady state if requested
1903	if(ice_update)then
1904
1905	do ig = 1, ngrid
1906
1907	delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
1908
1909	! add one hundred years/timesteps of evolution
1910	qsurf_hist(ig,igcm_h2o_ice) = &
1911	qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0
1912
1913	! if ice has gone -ve, set to zero
1914	if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
1915	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1916	qsurf_hist(ig,igcm_h2o_vap) = 0.0
1917	endif
1918
1919	! set maximum to ice
1920	if(qsurf_hist(ig,igcm_h2o_ice).gt.1.e3)then
1921	qsurf_hist(ig,igcm_h2o_ice) = 1.e3
1922	qsurf_hist(ig,igcm_h2o_vap) = 0.0
1923	endif
1924
1925	enddo
1926
1927	! enforce ice conservation
1928	ice_tot=0.0
1929	do ig = 1, ngrid
1930	ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
1931	enddo
1932	do ig = 1, ngrid
1933	qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
1934	enddo
1935
1936
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, &
1943	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1944	area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1945	cloudfrac,totcloudfrac,hice)
1946	endif
1947
1948	! -----------------------------------------------------------------
1949	! Saving statistics :
1950	! -----------------------------------------------------------------
1951	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1952	! which can later be used to make the statistic files of the run:
1953	! "stats") only possible in 3D runs !
1954
1955
1956	if (callstats) then
1957
1958	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1959	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1960	call wstats(ngrid,"fluxsurf_lw", &
1961	"Thermal IR radiative flux to surface","W.m-2",2, &
1962	fluxsurf_lw)
1963	! call wstats(ngrid,"fluxsurf_sw", &
1964	! "Solar radiative flux to surface","W.m-2",2, &
1965	! fluxsurf_sw_tot)
1966	call wstats(ngrid,"fluxtop_lw", &
1967	"Thermal IR radiative flux to space","W.m-2",2, &
1968	fluxtop_lw)
1969	! call wstats(ngrid,"fluxtop_sw", &
1970	! "Solar radiative flux to space","W.m-2",2, &
1971	! fluxtop_sw_tot)
1972	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1973	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1974	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1975	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1976	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1977
1978	if (tracer) then
1979	if (water) then
1980	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
1981	call wstats(ngrid,"vmr_h2ovapor", &
1982	"H2O vapour volume mixing ratio","mol/mol", &
1983	3,vmr)
1984	endif ! of if (water)
1985
1986	endif !tracer
1987
1988	if(lastcall) then
1989	write (,) "Writing stats..."
1990	call mkstats(ierr)
1991	endif
1992	endif !if callstats
1993
1994
1995	! ----------------------------------------------------------------------
1996	! output in netcdf file "DIAGFI", containing any variable for diagnostic
1997	! (output with period "ecritphy", set in "run.def")
1998	! ----------------------------------------------------------------------
1999	! writediagfi can also be called from any other subroutine for any variable.
2000	! but its preferable to keep all the calls in one place...
2001
2002	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2003	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2004
2005	if(1.eq.2)then
2006	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2007	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2008	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2009	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2010	call writediagfi(ngrid,'p','Pressure','Pa',3,pplay)
2011	endif
2012
2013	! Total energy balance diagnostics
2014	if(callrad.and.(.not.newtonian))then
2015	call writediagfi(ngrid,'ALB','Surface albedo',' ',2,albedo)
2016	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2017	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2018	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2019	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2020	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2021	!call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
2022	!call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
2023	!call writediagfi(ngrid,"vdifcdE","heat from vdifc surface","W m-2",2,vdifcdE)
2024	endif
2025
2026	! Temporary inclusions for heating diagnostics
2027	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2028	! call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2029	! call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2030	! call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2031
2032	! debugging
2033	!call writediagfi(ngrid,"vdifNC","H2O loss vdifc","kg m-2 s-1",2,vdifcncons)
2034	!call writediagfi(ngrid,"cadjNC","H2O loss convadj","kg m-2 s-1",2,cadjncons)
2035	call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2036	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2037	!print*,vdifcncons
2038	!call abort
2039
2040	! Output aerosols
2041	! call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
2042	! call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
2043	! call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
2044	! call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
2045
2046	! Output tracers
2047	if (tracer) then
2048	do iq=1,nq
2049	! call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2050	! call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2051	! 'kg m^-2',2,qsurf(1,iq) )
2052	call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2053	'kg m^-2',2,qsurf_hist(1,iq) )
2054	call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2055	'kg m^-2',2,qcol(1,iq) )
2056
2057	if(water)then
2058	call writediagfi(ngridmx,"H2Omaxcol","max. poss. H2O column","kg m^-2",2,H2Omaxcol)
2059	endif
2060
2061	if(watercond)then
2062	!call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2063	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2064	endif
2065
2066	if(waterrain)then
2067	call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
2068	call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
2069	endif
2070
2071	if(hydrology)then
2072	call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
2073	endif
2074
2075	if(ice_update)then
2076	call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
2077	call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
2078	endif
2079
2080	do ig=1,ngrid
2081	if(tau_col(ig).gt.1.e3)then
2082	print*,'WARNING: tau_col=',tau_col(ig)
2083	print*,'at ig=',ig,'in PHYSIQ'
2084	!call abort
2085	endif
2086	end do
2087
2088	call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2089
2090	enddo
2091	endif
2092
2093	! ----------------------------------------------------------
2094	! Output in netcdf file "diagsoil.nc" for subterranean
2095	! variables (output every "ecritphy", as for writediagfi)
2096	! ----------------------------------------------------------
2097
2098	! Write soil temperature
2099	!call writediagsoil(ngrid,"soiltemp","Soil temperature","K",3,tsoil)
2100
2101	else ! if(ngrid.eq.1)
2102
2103	! ----------------------------------------------------------------------
2104	! Output in grads file "g1d" (ONLY when using testphys1d)
2105	! ----------------------------------------------------------------------
2106
2107	if(countG1D.eq.saveG1D)then
2108
2109	call writeg1d(ngrid,1,tsurf,'tsurf','K')
2110	call writeg1d(ngrid,1,ps,'ps','Pa')
2111	call writeg1d(ngrid,nlayer,zt,'T','K')
2112	call writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg')
2113	call writeg1d(ngrid,nlayer,aerosol,'aerosol','SI')
2114	call writeg1d(ngrid,nlayer,reffrad,'reffrad','SI')
2115	call writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI')
2116	call writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI')
2117	call writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI')
2118
2119	! radiation balance
2120	call writeg1d(ngrid,1,ISR/totarea,'ISR','W m-2')
2121	call writeg1d(ngrid,1,ASR/totarea,'ASR','W m-2')
2122	call writeg1d(ngrid,1,OLR/totarea,'OLR','W m-2')
2123
2124	if(tracer) then
2125	do iq=1,nq
2126	call writeg1d(ngrid,1,qsurf(1,iq),trim(noms(iq))//'_s','kg m^-2')
2127	call writeg1d(ngrid,1,qcol(1,iq),trim(noms(iq))//'_c','kg m^-2')
2128	call writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
2129	end do
2130
2131	if(waterrain)then
2132	CALL writeg1d(ngrid,1,zdqsrain,'rainfall','kg m-2 s-1')
2133	CALL writeg1d(ngrid,1,zdqssnow,'snowfall','kg m-2 s-1')
2134	endif
2135
2136	end if
2137
2138	countG1D=1
2139	else
2140	countG1D=countG1D+1
2141	endif ! if time to save
2142
2143	endif ! if(ngrid.ne.1)
2144
2145	icount=icount+1
2146
2147	return
2148	end subroutine physiq

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

Download in other formats:

Original Format