Context Navigation

physiq.F90 @ 735

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

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

Download in other formats:

Original Format