Context Navigation

physiq.F90 @ 764

Last change on this file since 764 was 751, checked in by aslmd, 13 years ago
GENERIC. Sorry. Erased commits 749 and 750. Waiting for collegial approval. Only made a minor correction: reffrad1 is not useful in physiq and might cause problems with some compilers so it is removed
Property svn:executable set to ``*
File size: 70.5 KB

Line
1	subroutine physiq(ngrid,nlayer,nq, &
2	firstcall,lastcall, &
3	pday,ptime,ptimestep, &
4	pplev,pplay,pphi, &
5	pu,pv,pt,pq, &
6	pw, &
7	pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
8
9	use radinc_h, only : 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
647	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
648	ztim1,ztim2,ztim3,mu0,fract)
649
650	elseif (diurnal) then
651	ztim1=SIN(declin)
652	ztim2=COS(declin)COS(2.pi*(zday-.5))
653	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
654
655	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
656	ztim1,ztim2,ztim3,mu0,fract)
657
658	else
659
660	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
661	! WARNING: this function appears not to work in 1D
662
663	endif
664
665	if (corrk) then
666
667	! a) Call correlated-k radiative transfer scheme
668	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
669
670	if(kastprof)then
671	print*,'kastprof should not = true here'
672	call abort
673	endif
674	muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
675
676	! standard callcorrk
677	clearsky=.false.
678	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
679	albedo,emis,mu0,pplev,pplay,pt, &
680	tsurf,fract,dist_star,aerosol,muvar, &
681	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
682	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
683	reffrad,tau_col,cloudfrac,totcloudfrac, &
684	clearsky,firstcall,lastcall)
685
686	! Option to call scheme once more for clear regions
687	if(CLFvarying)then
688
689	! ---> PROBLEMS WITH ALLOCATED ARRAYS
690	! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
691	clearsky=.true.
692	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
693	albedo,emis,mu0,pplev,pplay,pt, &
694	tsurf,fract,dist_star,aerosol,muvar, &
695	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
696	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, &
697	reffrad,tau_col1,cloudfrac,totcloudfrac, &
698	clearsky,firstcall,lastcall)
699	clearsky = .false. ! just in case.
700
701	! Sum the fluxes and heating rates from cloudy/clear cases
702	do ig=1,ngrid
703	tf=totcloudfrac(ig)
704	ntf=1.-tf
705
706	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
707	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
708	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
709	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
710	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
711
712	zdtlw(ig,1:nlayermx) = ntfzdtlw1(ig,1:nlayermx) + tfzdtlw(ig,1:nlayermx)
713	zdtsw(ig,1:nlayermx) = ntfzdtsw1(ig,1:nlayermx) + tfzdtsw(ig,1:nlayermx)
714
715	OSR_nu(ig,1:L_NSPECTV) = ntfOSR_nu1(ig,1:L_NSPECTV) + tfOSR_nu(ig,1:L_NSPECTV)
716	OLR_nu(ig,1:L_NSPECTI) = ntfOLR_nu1(ig,1:L_NSPECTI) + tfOLR_nu(ig,1:L_NSPECTI)
717
718	enddo
719
720	endif !CLFvarying
721
722	! Radiative flux from the sky absorbed by the surface (W.m-2)
723	GSR=0.0
724	fluxrad_sky(1:ngridmx)=emis(1:ngridmx)fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)(1.-albedo(1:ngridmx))
725
726	if(noradsurf)then ! no lower surface; SW flux just disappears
727	GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
728	fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
729	print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
730	endif
731
732	! Net atmospheric radiative heating rate (K.s-1)
733	dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
734
735	elseif(newtonian)then
736
737	! b) Call Newtonian cooling scheme
738	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
739	call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
740
741	zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
742	! e.g. surface becomes proxy for 1st atmospheric layer ?
743
744	else
745
746	! c) Atmosphere has no radiative effect
747	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
748	fluxtop_dn(1:ngridmx) = fract(1:ngridmx)mu0(1:ngridmx)Fat1AU/dist_star**2
749	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
750	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
751	endif
752	fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
753	fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
754	fluxtop_lw(1:ngridmx) = emis(1:ngridmx)sigmatsurf(1:ngridmx)**4
755	! radiation skips the atmosphere entirely
756
757
758	dtrad(1:ngridmx,1:nlayermx)=0.0
759	! hence no atmospheric radiative heating
760
761	endif ! if corrk
762
763	endif ! of if(mod(icount-1,iradia).eq.0)
764
765
766	! Transformation of the radiative tendencies
767	! ------------------------------------------
768
769	zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
770	zplanck(1:ngridmx)=emis(1:ngridmx)sigmazplanck(1:ngridmx)*zplanck(1:ngridmx)
771	fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
772	pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
773
774	!-------------------------
775	! test energy conservation
776	if(enertest)then
777	dEtotSW = cppSUM(massarea(:,:)zdtsw(:,:))/totarea
778	dEtotLW = cppSUM(massarea(:,:)zdtlw(:,:))/totarea
779	dEtotsSW = SUM(fluxsurf_sw(:)(1.-albedo(:))area(:))/totarea
780	dEtotsLW = SUM((fluxsurf_lw(:)emis(:)-zplanck(:))area(:))/totarea
781	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
782	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
783	print*,'---------------------------------------------------------------'
784	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
785	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
786	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
787	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
788	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
789	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
790	endif
791	!-------------------------
792
793	endif ! of if (callrad)
794
795	!-----------------------------------------------------------------------
796	! 4. Vertical diffusion (turbulent mixing):
797	! -----------------------------------------
798
799	if (calldifv) then
800
801	zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
802
803	zdum1(1:ngridmx,1:nlayermx)=0.0
804	zdum2(1:ngridmx,1:nlayermx)=0.0
805
806
807	!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff
808	if (UseTurbDiff) then
809
810	call turbdiff(ngrid,nlayer,nq,rnat, &
811	ptimestep,capcal,lwrite, &
812	pplay,pplev,zzlay,zzlev,z0, &
813	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
814	zdum1,zdum2,pdt,pdq,zflubid, &
815	zdudif,zdvdif,zdtdif,zdtsdif, &
816	sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
817
818	else
819
820	zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
821
822	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
823	ptimestep,capcal,lwrite, &
824	pplay,pplev,zzlay,zzlev,z0, &
825	pu,pv,zh,pq,tsurf,emis,qsurf, &
826	zdum1,zdum2,zdh,pdq,zflubid, &
827	zdudif,zdvdif,zdhdif,zdtsdif, &
828	sensibFlux,q2,zdqdif,zdqsdif,lastcall)
829
830	zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
831	zdqevap(1:ngridmx,1:nlayermx)=0.
832
833	end if !turbdiff
834
835	pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
836	pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
837	pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
838	zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
839	if (tracer) then
840	pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx)
841	dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
842	end if ! of if (tracer)
843
844	!-------------------------
845	! test energy conservation
846	if(enertest)then
847	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
848	do ig = 1, ngrid
849	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
850	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
851	enddo
852	dEtot = SUM(dEdiff(:)*area(:))/totarea
853	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
854	dEtots = SUM(dEdiffs(:)*area(:))/totarea
855	AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
856	if (UseTurbDiff) then
857	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
858	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
859	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
860	else
861	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
862	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
863	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
864	end if
865	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
866	! but not given back elsewhere
867	endif
868	!-------------------------
869
870	!-------------------------
871	! test water conservation
872	if(watertest.and.water)then
873	dWtot = SUM(massarea(:,:)zdqdif(:,:,igcm_h2o_vap))ptimestep/totarea
874	dWtots = SUM(zdqsdif(:,igcm_h2o_vap)area(:))ptimestep/totarea
875	do ig = 1, ngrid
876	vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
877	Enddo
878	dWtot = dWtot + SUM(massarea(:,:)zdqdif(:,:,igcm_h2o_ice))ptimestep/totarea
879	dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)area(:))ptimestep/totarea
880	do ig = 1, ngrid
881	vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
882	Enddo
883	nconsMAX=MAXVAL(vdifcncons(:))
884
885	print*,'---------------------------------------------------------------'
886	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
887	print*,'In difv surface water change =',dWtots,' kg m-2'
888	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
889	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
890
891	endif
892	!-------------------------
893
894	else
895
896	if(.not.newtonian)then
897
898	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
899
900	endif
901
902	endif ! of if (calldifv)
903
904
905	!-----------------------------------------------------------------------
906	! 5. Dry convective adjustment:
907	! -----------------------------
908
909	if(calladj) then
910
911	zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
912	zduadj(1:ngridmx,1:nlayermx)=0.0
913	zdvadj(1:ngridmx,1:nlayermx)=0.0
914	zdhadj(1:ngridmx,1:nlayermx)=0.0
915
916
917	call convadj(ngrid,nlayer,nq,ptimestep, &
918	pplay,pplev,zpopsk, &
919	pu,pv,zh,pq, &
920	pdu,pdv,zdh,pdq, &
921	zduadj,zdvadj,zdhadj, &
922	zdqadj)
923
924	pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
925	pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
926	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
927	zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
928
929	if(tracer) then
930	pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx)
931	end if
932
933	!-------------------------
934	! test energy conservation
935	if(enertest)then
936	dEtot=cppSUM(massarea(:,:)zdtadj(:,:))/totarea
937	print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
938	endif
939	!-------------------------
940
941	!-------------------------
942	! test water conservation
943	if(watertest)then
944	dWtot = SUM(massarea(:,:)zdqadj(:,:,igcm_h2o_vap))ptimestep/totarea
945	do ig = 1, ngrid
946	cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
947	Enddo
948	dWtot = dWtot + SUM(massarea(:,:)zdqadj(:,:,igcm_h2o_ice))ptimestep/totarea
949	do ig = 1, ngrid
950	cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
951	Enddo
952	nconsMAX=MAXVAL(cadjncons(:))
953
954	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
955	print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
956	endif
957	!-------------------------
958
959	endif ! of if(calladj)
960
961	!-----------------------------------------------------------------------
962	! 6. Carbon dioxide condensation-sublimation:
963	! -------------------------------------------
964
965	if (co2cond) then
966	if (.not.tracer) then
967	print*,'We need a CO2 ice tracer to condense CO2'
968	call abort
969	endif
970	print*, 'we are in co2cond!!!!!'
971	call condense_cloud(ngrid,nlayer,nq,ptimestep, &
972	capcal,pplay,pplev,tsurf,pt, &
973	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
974	qsurf(1,igcm_co2_ice),albedo,emis, &
975	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
976	zdqc,reffrad)
977
978
979	pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
980	pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
981	pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
982	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
983
984	pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx)
985	! Note: we do not add surface co2ice tendency
986	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
987
988	!-------------------------
989	! test energy conservation
990	if(enertest)then
991	dEtot = cppSUM(massarea(:,:)zdtc(:,:))/totarea
992	dEtots = SUM(capcal(:)zdtsurfc(:)area(:))/totarea
993	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
994	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
995	endif
996	!-------------------------
997
998	endif ! of if (co2cond)
999
1000
1001	!-----------------------------------------------------------------------
1002	! 7. Specific parameterizations for tracers
1003	! -----------------------------------------
1004
1005	if (tracer) then
1006
1007	! 7a. Water and ice
1008	! ---------------
1009	if (water) then
1010
1011	! ----------------------------------------
1012	! Water ice condensation in the atmosphere
1013	! ----------------------------------------
1014	if(watercond.and.(RLVTT.gt.1.e-8))then
1015
1016	! ----------------
1017	! Moist convection
1018	! ----------------
1019	! Re-evaporate cloud water/ice
1020	! dqevap1(1:ngridmx,1:nlayermx)=0.
1021	! dtevap1(1:ngridmx,1:nlayermx)=0.
1022	! call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
1023
1024	dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
1025	dtmoist(1:ngridmx,1:nlayermx)=0.
1026
1027	call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1028
1029	pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
1030	+dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1031	pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
1032	+dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1033	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
1034
1035	! do l=1,nlayermx
1036	! do ig=1,ngridmx
1037	! if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
1038	! if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
1039	! if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l
1040	! if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
1041	! enddo
1042	! enddo
1043	! print,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)86400.
1044
1045	!-------------------------
1046	! test energy conservation
1047	if(enertest)then
1048	dEtot=cppSUM(massarea(:,:)dtmoist(:,:))/totarea
1049	do ig = 1, ngrid
1050	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
1051	enddo
1052	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
1053	print,'In moistadj MAX atmospheric energy change =',MAXVAL(pdt(:,:))86400.,'K/day'
1054
1055	! test energy conservation
1056	dWtot = SUM(massarea(:,:)dqmoist(:,:,igcm_h2o_vap))ptimestep/totarea
1057	dWtot = dWtot + SUM(massarea(:,:)dqmoist(:,:,igcm_h2o_ice))ptimestep/totarea
1058	print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
1059	endif
1060	!-------------------------
1061
1062
1063	! --------------------------------
1064	! Large scale condensation/evaporation
1065	! --------------------------------
1066
1067	call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
1068
1069	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
1070	pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
1071	pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
1072
1073	!-------------------------
1074	! test energy conservation
1075	if(enertest)then
1076	do ig = 1, ngrid
1077	lscaledE(ig) = cppSUM(mass(:,:)dtlscale(:,:))
1078	enddo
1079	dEtot=cppSUM(massarea(:,:)(dtlscale(:,:)))/totarea
1080	if(isnan(dEtot)) then
1081	print*,'Nan in largescale, abort'
1082	STOP
1083	endif
1084	print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
1085
1086	! test water conservation
1087	dWtot = SUM(massarea(:,:)dqvaplscale(:,:))ptimestep/totarea
1088	dWtot = dWtot + SUM(massarea(:,:)dqcldlscale(:,:))ptimestep/totarea
1089	print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
1090	endif
1091	!-------------------------
1092
1093	! compute cloud fraction
1094	do l = 1, nlayer
1095	do ig = 1,ngrid
1096	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1097	enddo
1098	enddo
1099
1100	endif ! of if (watercondense)
1101
1102
1103	! --------------------------------
1104	! Water ice / liquid precipitation
1105	! --------------------------------
1106	if(waterrain)then
1107
1108	zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
1109	zdqsrain(1:ngridmx) = 0.0
1110	zdqssnow(1:ngridmx) = 0.0
1111
1112	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
1113	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
1114
1115	pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
1116	+zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1117	pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
1118	+zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1119	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
1120	dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
1121	dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
1122	rainout(1:ngridmx) = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic
1123
1124
1125	!-------------------------
1126	! test energy conservation
1127	if(enertest)then
1128	dEtot=cppSUM(massarea(:,:)zdtrain(:,:))/totarea
1129	print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
1130	dItot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_ice))/totareaRLVTT/cpp
1131	dItot = dItot + SUM(area(:)zdqssnow(:))/totareaRLVTT/cpp
1132	dVtot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_vap))ptimestep/totarea
1133	dVtot = dVtot + SUM(area(:)zdqsrain(:))/totareaRLVTT/cpp
1134	dEtot = dItot + dVtot
1135	print*,'In rain dItot =',dItot,' W m-2'
1136	print*,'In rain dVtot =',dVtot,' W m-2'
1137	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
1138
1139	! test water conservation
1140	dWtot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_vap))ptimestep/totarea
1141	dWtot = dWtot + SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_ice))ptimestep/totarea
1142	dWtots = SUM((zdqsrain(:)+zdqssnow(:))area(:))ptimestep/totarea
1143	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
1144	print*,'In rain surface water change =',dWtots,' kg m-2'
1145	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
1146	endif
1147	!-------------------------
1148
1149	end if ! of if (waterrain)
1150	end if ! of if (water)
1151
1152
1153	! 7c. Aerosol particles
1154	! -------------------
1155	! -------------
1156	! Sedimentation
1157	! -------------
1158	if (sedimentation) then
1159	zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
1160	zdqssed(1:ngridmx,1:nqmx) = 0.0
1161
1162
1163	!-------------------------
1164	! find qtot
1165	if(watertest)then
1166	iq=3
1167	dWtot = SUM(massarea(:,:)pq(:,:,iq))ptimestep/totarea
1168	dWtots = SUM(massarea(:,:)pdq(:,:,iq))ptimestep/totarea
1169	print*,'Before sedim pq =',dWtot,' kg m-2'
1170	print*,'Before sedim pdq =',dWtots,' kg m-2'
1171	endif
1172	!-------------------------
1173
1174	call callsedim(ngrid,nlayer,ptimestep, &
1175	pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
1176
1177	!-------------------------
1178	! find qtot
1179	if(watertest)then
1180	iq=3
1181	dWtot = SUM(massarea(:,:)pq(:,:,iq))ptimestep/totarea
1182	dWtots = SUM(massarea(:,:)pdq(:,:,iq))ptimestep/totarea
1183	print*,'After sedim pq =',dWtot,' kg m-2'
1184	print*,'After sedim pdq =',dWtots,' kg m-2'
1185	endif
1186	!-------------------------
1187
1188	! for now, we only allow H2O ice to sediment
1189	! and as in rain.F90, whether it falls as rain or snow depends
1190	! only on the surface temperature
1191	pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
1192	dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
1193
1194	!-------------------------
1195	! test water conservation
1196	if(watertest)then
1197	dWtot = SUM(massarea(:,:)(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))ptimestep/totarea
1198	dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))area(:))ptimestep/totarea
1199	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
1200	print*,'In sedim surface ice change =',dWtots,' kg m-2'
1201	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
1202	endif
1203	!-------------------------
1204
1205	end if ! of if (sedimentation)
1206
1207
1208	! 7d. Updates
1209	! ---------
1210
1211	! -----------------------------------
1212	! Updating atm mass and tracer budget
1213	! -----------------------------------
1214
1215	if(mass_redistrib) then
1216
1217	zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
1218	( zdqevap(1:ngridmx,1:nlayermx) &
1219	+ zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
1220	+ dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
1221	+ dqvaplscale(1:ngridmx,1:nlayermx) )
1222
1223
1224	call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
1225	call writediagfi(ngridmx,"mass","mass"," ",3,mass)
1226
1227	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1228	rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1229	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1230	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1231
1232
1233	pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
1234	dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
1235	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
1236	pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
1237	pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
1238	pdpsrf(1:ngridmx) = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
1239	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
1240
1241	! print,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimesteppdq(211,1:nlayermx,igcm_h2o_vap)
1242	endif
1243
1244
1245
1246	! ---------------------------------
1247	! Updating tracer budget on surface
1248	! ---------------------------------
1249
1250	if(hydrology)then
1251
1252	call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
1253	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
1254	! note: for now, also changes albedo in the subroutine
1255
1256	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
1257	qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
1258	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
1259
1260	!-------------------------
1261	! test energy conservation
1262	if(enertest)then
1263	dEtots = SUM(area(:)capcal(:)zdtsurf_hyd(:))/totarea
1264	print*,'In hydrol surface energy change =',dEtots,' W m-2'
1265	endif
1266	!-------------------------
1267
1268	!-------------------------
1269	! test water conservation
1270	if(watertest)then
1271	dWtots = SUM(dqs_hyd(:,igcm_h2o_ice)area(:))ptimestep/totarea
1272	print*,'In hydrol surface ice change =',dWtots,' kg m-2'
1273	dWtots = SUM(dqs_hyd(:,igcm_h2o_vap)area(:))ptimestep/totarea
1274	print*,'In hydrol surface water change =',dWtots,' kg m-2'
1275	print*,'---------------------------------------------------------------'
1276	endif
1277	!-------------------------
1278
1279	ELSE ! of if (hydrology)
1280
1281	qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
1282
1283	END IF ! of if (hydrology)
1284
1285	! Add qsurf to qsurf_hist, which is what we save in
1286	! diagfi.nc etc. At the same time, we set the water
1287	! content of ocean gridpoints back to zero, in order
1288	! to avoid rounding errors in vdifc, rain
1289	qsurf_hist(:,:) = qsurf(:,:)
1290
1291	if(ice_update)then
1292	ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
1293	endif
1294
1295	endif ! of if (tracer)
1296
1297	!-----------------------------------------------------------------------
1298	! 9. Surface and sub-surface soil temperature
1299	!-----------------------------------------------------------------------
1300
1301
1302	! Increment surface temperature
1303	tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx)
1304
1305	! Compute soil temperatures and subsurface heat flux
1306	if (callsoil) then
1307	call soil(ngrid,nsoilmx,.false.,inertiedat, &
1308	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1309	endif
1310
1311	!-------------------------
1312	! test energy conservation
1313	if(enertest)then
1314	dEtots = SUM(area(:)capcal(:)zdtsurf(:))/totarea
1315	print*,'Surface energy change =',dEtots,' W m-2'
1316	endif
1317	!-------------------------
1318
1319	!-----------------------------------------------------------------------
1320	! 10. Perform diagnostics and write output files
1321	!-----------------------------------------------------------------------
1322
1323	! -------------------------------
1324	! Dynamical fields incrementation
1325	! -------------------------------
1326	! For output only: the actual model integration is performed in the dynamics
1327
1328	! temperature, zonal and meridional wind
1329	zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
1330	zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
1331	zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
1332
1333	! diagnostic
1334	zdtdyn(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
1335	ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
1336
1337	if(firstcall)then
1338	zdtdyn(1:ngridmx,1:nlayermx)=0.0
1339	endif
1340
1341	! dynamical heating diagnostic
1342	do ig=1,ngrid
1343	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp/ptimestep
1344	enddo
1345
1346	! tracers
1347	zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
1348
1349	! surface pressure
1350	ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
1351
1352	! pressure
1353	do l=1,nlayer
1354	zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
1355	zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
1356	enddo
1357
1358	! ---------------------------------------------------------
1359	! Surface and soil temperature information
1360	! ---------------------------------------------------------
1361
1362	Ts1 = SUM(area(:)*tsurf(:))/totarea
1363	Ts2 = MINVAL(tsurf(:))
1364	Ts3 = MAXVAL(tsurf(:))
1365	if(callsoil)then
1366	TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1367	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1368	print*,Ts1,Ts2,Ts3,TsS
1369	else
1370	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1371	print*,Ts1,Ts2,Ts3
1372	endif
1373
1374	! ---------------------------------------------------------
1375	! Check the energy balance of the simulation during the run
1376	! ---------------------------------------------------------
1377
1378	if(corrk)then
1379
1380	ISR = SUM(area(:)*fluxtop_dn(:))/totarea
1381	ASR = SUM(area(:)*fluxabs_sw(:))/totarea
1382	OLR = SUM(area(:)*fluxtop_lw(:))/totarea
1383	GND = SUM(area(:)*fluxgrd(:))/totarea
1384	DYN = SUM(area(:)*fluxdyn(:))/totarea
1385	do ig=1,ngrid
1386	if(fluxtop_dn(ig).lt.0.0)then
1387	print*,'fluxtop_dn has gone crazy'
1388	print*,'fluxtop_dn=',fluxtop_dn(ig)
1389	print*,'tau_col=',tau_col(ig)
1390	print*,'aerosol=',aerosol(ig,:,:)
1391	print*,'temp= ',pt(ig,:)
1392	print*,'pplay= ',pplay(ig,:)
1393	call abort
1394	endif
1395	end do
1396
1397	if(ngridmx.eq.1)then
1398	DYN=0.0
1399	endif
1400
1401	print*,' ISR ASR OLR GND DYN [W m^-2]'
1402	print*, ISR,ASR,OLR,GND,DYN
1403
1404	if(enertest)then
1405	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1406	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1407	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1408	endif
1409
1410	if(meanOLR)then
1411	if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
1412	! to record global radiative balance
1413	open(92,file="rad_bal.out",form='formatted',position='append')
1414	write(92,*) zday,ISR,ASR,OLR
1415	close(92)
1416	open(93,file="tem_bal.out",form='formatted',position='append')
1417	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1418	close(93)
1419	endif
1420	endif
1421
1422	endif
1423
1424
1425	! ------------------------------------------------------------------
1426	! Diagnostic to test radiative-convective timescales in code
1427	! ------------------------------------------------------------------
1428	if(testradtimes)then
1429	open(38,file="tau_phys.out",form='formatted',position='append')
1430	ig=1
1431	do l=1,nlayer
1432	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1433	enddo
1434	close(38)
1435	print*,'As testradtimes enabled,'
1436	print*,'exiting physics on first call'
1437	call abort
1438	endif
1439
1440	! ---------------------------------------------------------
1441	! Compute column amounts (kg m-2) if tracers are enabled
1442	! ---------------------------------------------------------
1443	if(tracer)then
1444	qcol(1:ngridmx,1:nqmx)=0.0
1445	do iq=1,nq
1446	do ig=1,ngridmx
1447	qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
1448	enddo
1449	enddo
1450
1451	! Generalised for arbitrary aerosols now. (LK)
1452	reffcol(1:ngridmx,1:naerkind)=0.0
1453	if(co2cond.and.(iaero_co2.ne.0))then
1454	call co2_reffrad(zq,reffrad)
1455	do ig=1,ngridmx
1456	reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)reffrad(ig,1:nlayermx,iaero_co2)mass(ig,1:nlayermx))
1457	enddo
1458	endif
1459	if(water.and.(iaero_h2o.ne.0))then
1460	call h2o_reffrad(zq,zt,reffrad,nueffrad)
1461	do ig=1,ngridmx
1462	reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)reffrad(ig,1:nlayermx,iaero_h2o)mass(ig,1:nlayermx))
1463	enddo
1464	endif
1465
1466	endif
1467
1468	! ---------------------------------------------------------
1469	! Test for water conservation if water is enabled
1470	! ---------------------------------------------------------
1471
1472	if(water)then
1473
1474	icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
1475	liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
1476	icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
1477	vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
1478
1479	h2otot = icesrf + liqsrf + icecol + vapcol
1480
1481	print*,' Total water amount [kg m^-2]: ',h2otot
1482	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
1483	print*, icesrf,liqsrf,icecol,vapcol
1484
1485	if(meanOLR)then
1486	if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
1487	! to record global water balance
1488	open(98,file="h2o_bal.out",form='formatted',position='append')
1489	write(98,*) zday,icesrf,liqsrf,icecol,vapcol
1490	close(98)
1491	endif
1492	endif
1493
1494	endif
1495
1496	! ---------------------------------------------------------
1497	! Calculate RH for diagnostic if water is enabled
1498	! ---------------------------------------------------------
1499
1500	if(water)then
1501	do l = 1, nlayer
1502	do ig = 1, ngrid
1503	! call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
1504	call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
1505
1506	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
1507	enddo
1508	enddo
1509
1510	! compute maximum possible H2O column amount (100% saturation)
1511	do ig=1,ngrid
1512	H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
1513	enddo
1514
1515	endif
1516
1517
1518	print*,''
1519	print,'--> Ls =',zls180./pi
1520	! -------------------------------------------------------------------
1521	! Writing NetCDF file "RESTARTFI" at the end of the run
1522	! -------------------------------------------------------------------
1523	! Note: 'restartfi' is stored just before dynamics are stored
1524	! in 'restart'. Between now and the writting of 'restart',
1525	! there will have been the itau=itau+1 instruction and
1526	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1527	! thus we store for time=time+dtvr
1528
1529	if(lastcall) then
1530	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1531
1532
1533	! Update surface ice distribution to iterate to steady state if requested
1534	if(ice_update)then
1535
1536	do ig = 1, ngrid
1537
1538	delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
1539
1540	! add multiple years of evolution
1541	qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
1542
1543	! if ice has gone -ve, set to zero
1544	if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
1545	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1546	endif
1547
1548	! if ice is seasonal, set to zero (NEW)
1549	if(ice_min(ig).lt.0.01)then
1550	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1551	endif
1552
1553	enddo
1554
1555	! enforce ice conservation
1556	ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
1557	qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
1558
1559	endif
1560
1561	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1562	call physdem1("restartfi.nc",long,lati,nsoilmx,nq, &
1563	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1564	area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1565	cloudfrac,totcloudfrac,hice)
1566	endif
1567
1568	! -----------------------------------------------------------------
1569	! Saving statistics :
1570	! -----------------------------------------------------------------
1571	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1572	! which can later be used to make the statistic files of the run:
1573	! "stats") only possible in 3D runs !
1574
1575
1576	if (callstats) then
1577
1578	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1579	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1580	call wstats(ngrid,"fluxsurf_lw", &
1581	"Thermal IR radiative flux to surface","W.m-2",2, &
1582	fluxsurf_lw)
1583	! call wstats(ngrid,"fluxsurf_sw", &
1584	! "Solar radiative flux to surface","W.m-2",2, &
1585	! fluxsurf_sw_tot)
1586	call wstats(ngrid,"fluxtop_lw", &
1587	"Thermal IR radiative flux to space","W.m-2",2, &
1588	fluxtop_lw)
1589	! call wstats(ngrid,"fluxtop_sw", &
1590	! "Solar radiative flux to space","W.m-2",2, &
1591	! fluxtop_sw_tot)
1592
1593	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1594	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1595	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1596
1597	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1598	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1599	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1600	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1601	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1602
1603	if (tracer) then
1604	do iq=1,nq
1605	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1606	call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1607	'kg m^-2',2,qsurf(1,iq) )
1608
1609	call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
1610	'kg m^-2',2,qcol(1,iq) )
1611	! call wstats(ngridmx,trim(noms(iq))//'_reff', &
1612	! trim(noms(iq))//'_reff', &
1613	! 'm',3,reffrad(1,1,iq))
1614	end do
1615	if (water) then
1616	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
1617	call wstats(ngrid,"vmr_h2ovapor", &
1618	"H2O vapour volume mixing ratio","mol/mol", &
1619	3,vmr)
1620	endif ! of if (water)
1621
1622	endif !tracer
1623
1624	if(lastcall) then
1625	write (,) "Writing stats..."
1626	call mkstats(ierr)
1627	endif
1628	endif !if callstats
1629
1630
1631	! ----------------------------------------------------------------------
1632	! output in netcdf file "DIAGFI", containing any variable for diagnostic
1633	! (output with period "ecritphy", set in "run.def")
1634	! ----------------------------------------------------------------------
1635	! writediagfi can also be called from any other subroutine for any variable.
1636	! but its preferable to keep all the calls in one place...
1637
1638	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1639	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1640	call writediagfi(ngrid,"temp","temperature","K",3,zt)
1641	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
1642	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1643	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1644	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1645	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
1646
1647	! Total energy balance diagnostics
1648	if(callrad.and.(.not.newtonian))then
1649	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
1650	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1651	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1652	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1653	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
1654	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1655	endif
1656
1657	if(enertest) then
1658	if (calldifv) then
1659	call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
1660	! call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
1661	! call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
1662	! call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
1663	call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
1664	endif
1665	if (corrk) then
1666	! call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
1667	! call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
1668	endif
1669	if(watercond) then
1670	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
1671	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
1672	call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
1673	! call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
1674	! call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
1675	endif
1676	endif
1677
1678	! Temporary inclusions for heating diagnostics
1679	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1680	! call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1681	! call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1682	! call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
1683
1684	! debugging
1685	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
1686	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
1687
1688	! Output aerosols
1689	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
1690	call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
1691	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
1692	call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
1693	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
1694	call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
1695	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
1696	call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
1697
1698	! Output tracers
1699	if (tracer) then
1700	do iq=1,nq
1701	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1702	! call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1703	! 'kg m^-2',2,qsurf(1,iq) )
1704	call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1705	'kg m^-2',2,qsurf_hist(1,iq) )
1706	call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
1707	'kg m^-2',2,qcol(1,iq) )
1708
1709	if(watercond.or.CLFvarying)then
1710	! call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
1711	! call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
1712	! call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
1713	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
1714	endif
1715
1716	if(waterrain)then
1717	call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
1718	call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
1719	endif
1720
1721	if(hydrology)then
1722	call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
1723	endif
1724
1725	if(ice_update)then
1726	call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
1727	call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
1728	endif
1729
1730	do ig=1,ngrid
1731	if(tau_col(ig).gt.1.e3)then
1732	print*,'WARNING: tau_col=',tau_col(ig)
1733	print*,'at ig=',ig,'in PHYSIQ'
1734	endif
1735	end do
1736
1737	call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
1738
1739	enddo
1740	endif
1741
1742	! output spectrum
1743	if(specOLR.and.corrk)then
1744	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
1745	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
1746	endif
1747
1748
1749	icount=icount+1
1750
1751	! deallocate gas variables
1752	if (lastcall) then
1753	IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
1754	IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac ) ! both allocated in su_gases.F90
1755	endif
1756
1757
1758	return
1759	end subroutine physiq

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

Download in other formats:

Original Format