Context Navigation

physiq.F90 @ 726

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

17/07/2012 == JL for LK

Generalization of aerosol scheme:
- any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
- addition of a module with the id numbers for aerosols (aerosol_mod.F90).
- initialization of aerosols id numbers in iniaerosol.F90
- compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

If starting from an old start file, recreate start file with the q=0 option in newstart.e.
update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

To add other aerosols, see Laura Kerber.

Property svn:executable set to *

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

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

Download in other formats:

Original Format