Context Navigation

physiq_mod.F90 @ 1644

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

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

Download in other formats:

Original Format