Context Navigation

physiq_mod.F90 @ 2101

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

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

Download in other formats:

Original Format