Context Navigation

physiq_mod.F90 @ 1834

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

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

Download in other formats:

Original Format