Context Navigation

physiq.F90 @ 1448

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

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

Download in other formats:

Original Format