Context Navigation

physiq_mod.F90 @ 3604

Last change on this file since 3604 was 3562, checked in by mmaurice, 2 months ago

Generic PCM

1D restart operational: a restart1D.nc file is created that contains
psurf, tracers, winds and temperature profiles. A retartfi.nc file is
also created. Move those to and start1D.nc and startfi.nc and set
"restart" flag to .true. in rcm1d.def to restart from the files (also
make sure that day0 corresponds to the value in startfi.nc).

Property svn:executable set to *

File size: 124.1 KB

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

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

Download in other formats:

Original Format