Context Navigation

physiq_mod.F90 @ 3083

Last change on this file since 3083 was 3081, checked in by yjaziri, 22 months ago

Generic PCM:

Move mu0 (cosinus of solar zenith angle) calculation out of callrad
This was a issue for photochemistry because mu0 was updated only iradia times

Property svn:executable set to *

File size: 118.4 KB

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

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

Download in other formats:

Original Format