Context Navigation

physiq_mod.F90 @ 2276

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

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

Download in other formats:

Original Format