Context Navigation

physiq_mod.F90 @ 1836

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

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

Download in other formats:

Original Format