Context Navigation

physiq_mod.F90 @ 1694

Last change on this file since 1694 was 1682, checked in by emillour, 8 years ago

All GCMs: set things up to enable pluging physics with dynamico

dyn3d
gcm.F90 : move I/O initialization (dates) to be done before physics

initialization

dyn3dpar
gcm.F : move I/O initialization (dates) to be done before physics

initialization

dynphy_lonlat:
inigeomphy_mod.F90 : add ind_cell_glo computation and transfer

to init_geometry

phy_common:
geometry_mod.F90 : add ind_cell_glo module variable to store global

column index

print_control_mod.F90 : make initialization occur via init_print_control_mod

to avoid circular module dependencies

init_print_control_mod.F90 : added to initialize print_control_mod module

variables

mod_phys_lmdz_mpi_data.F90 : use print_control_mod (rather than iniprint.h)
mod_phys_lmdz_para.F90 : use print_control_mod (rather than iniprint.h)
mod_phys_lmdz_omp_data.F90 : add is_omp_master (alias of is_omp_root) module

variable and use print_control_mod (rather than
iniprint.h)

physics_distribution_mod.F90 : add call to init_dimphy in

init_physics_distribution

xios_writefield.F90 : generic routine to output field with XIOS (for debug)

misc:
handle_err_m.F90 : call abort_physic, rather than abort_gcm
wxios.F90 : updates to enable unstructured grids

set module variable g_ctx_name to "LMDZ"
wxios_init(): remove call to wxios_context_init
wxios_context_init(): call xios_context_initialize with COMM_LMDZ_PHY
add routine wxios_set_context() to get handle and set context to XIOS
wxios_domain_param(): change arguments and generate the domain in-place
add wxios_domain_param_unstructured(): generate domain for unstructured case

NB: access is via "domain group" (whereas it is via "domain" in

wxios_domain_param)

dynphy_lonlat/phy[std|mars|venus|titan]:
iniphysiq_mod.F90 : Remove call to init_dimphy (which is now done in

phy_common/physics_distribution_mod.F90)

Property svn:executable set to *

File size: 94.1 KB

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

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

Download in other formats:

Original Format