Context Navigation

physiq.F90 @ 778

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

05/09/2012 == JL

Corrects typo in previous commit to account for possible eccentricity in

tlocked case

Property svn:executable set to *

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

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

Download in other formats:

Original Format