Context Navigation

physiq.F90 @ 1525

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

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

Download in other formats:

Original Format