Context Navigation

physiq_mod.F90 @ 2546

Last change on this file since 2546 was 2542, checked in by aslmd, 4 years ago

Generic GCM:

Large update of the chemical modules

Read chemical network from input files
Init chemistry with ModernTrac?
Photolysis online calculation

Property svn:executable set to *

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

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

Download in other formats:

Original Format