Context Navigation

physiq.F90 @ 1529

Last change on this file since 1529 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

Property svn:executable set to *

File size: 91.4 KB

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

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

Download in other formats:

Original Format