Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90 @ 224

Last change on this file since 224 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

Property svn:executable set to *

File size: 92.4 KB

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

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

Download in other formats: