Context Navigation

physiq_mod.F90 @ 2757

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

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

Download in other formats:

Original Format