Context Navigation

physiq_mod.F90 @ 2693

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

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

Download in other formats:

Original Format