Context Navigation

physiq_mod.F90 @ 2509

Last change on this file since 2509 was 2372, checked in by aslmd, 5 years ago

added a MESOSCALE flag at one location where it was missing (file_ok is not declared because also under the MESOSCALE pre-compiling flag)

Property svn:executable set to *
Property svn:mergeinfo set to
/trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 merged eligible

File size: 69.6 KB

Rev	Line
[1549]	1	module physiq_mod
	2
	3	implicit none
	4
	5	contains
	6
[751]	7	subroutine physiq(ngrid,nlayer,nq, &
[787]	8	nametrac, &
[253]	9	firstcall,lastcall, &
	10	pday,ptime,ptimestep, &
	11	pplev,pplay,pphi, &
	12	pu,pv,pt,pq, &
[1312]	13	flxw, &
[1576]	14	pdu,pdv,pdt,pdq,pdpsrf)
[253]	15
[1788]	16	use radinc_h, only : L_NSPECTI,L_NSPECTV
[2328]	17	use radcommon_h, only: sigma, gzlat, grav, BWNV
[1958]	18	use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
[1216]	19	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
[1327]	20	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
[1216]	21	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
[2242]	22	use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file
[1543]	23	use geometry_mod, only: latitude, longitude, cell_area
[1542]	24	USE comgeomfi_h, only: totarea, totarea_planet
[1795]	25	USE tracer_h
[1525]	26	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
[1670]	27	use phyetat0_mod, only: phyetat0
[1216]	28	use phyredem, only: physdem0, physdem1
[1295]	29	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
	30	use mod_phys_lmdz_para, only : is_master
[1308]	31	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
	32	obliquit, nres, z0
[1524]	33	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
	34	use time_phylmdz_mod, only: daysec
[2366]	35	#ifndef MESOSCALE
[1914]	36	use logic_mod, only: moyzon_ch
[1966]	37	use moyzon_mod, only: zphibar, zphisbar, zplevbar, zplaybar, &
	38	zzlevbar, zzlaybar, ztfibar, zqfibar
[2366]	39	#endif
[1397]	40	use callkeys_mod
[2291]	41	use phys_state_var_mod
	42	use turb_mod, only : q2,sensibFlux,turb_resolved
	43	#ifndef MESOSCALE
[1622]	44	use vertical_layers_mod, only: presnivs, pseudoalt
[1896]	45	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
[2291]	46	#else
	47	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
	48	comm_FLUXTOP_DN,comm_FLUXABS_SW,&
	49	comm_FLUXTOP_LW,comm_FLUXSURF_SW,&
	50	comm_FLUXSURF_LW,comm_FLXGRD
	51	#endif
[1623]	52	#ifdef CPP_XIOS
[1622]	53	use xios_output_mod, only: initialize_xios_output, &
	54	update_xios_timestep, &
	55	send_xios_field
[1896]	56	use wxios, only: wxios_context_init, xios_context_finalize
[1623]	57	#endif
[1926]	58	use muphy_diag
[253]	59	implicit none
	60
	61
	62	!==================================================================
	63	!
	64	! Purpose
	65	! -------
	66	! Central subroutine for all the physics parameterisations in the
	67	! universal model. Originally adapted from the Mars LMDZ model.
	68	!
	69	! The model can be run without or with tracer transport
	70	! depending on the value of "tracer" in file "callphys.def".
	71	!
	72	!
	73	! It includes:
	74	!
[1477]	75	! I. Initialization :
	76	! I.1 Firstcall initializations.
	77	! I.2 Initialization for every call to physiq.
[253]	78	!
[1477]	79	! II. Compute radiative transfer tendencies (longwave and shortwave) :
	80	! II.a Option 1 : Call correlated-k radiative transfer scheme.
	81	! II.b Option 2 : Call Newtonian cooling scheme.
	82	! II.c Option 3 : Atmosphere has no radiative effect.
	83	!
	84	! III. Vertical diffusion (turbulent mixing) :
	85	!
	86	! IV. Dry Convective adjustment :
	87	!
[1647]	88	! V. Tracers
[1926]	89	! V.1. Microphysics
	90	! V.2. Chemistry
[1672]	91	! V.3. Updates (pressure variations, surface budget).
	92	! V.4. Surface Tracer Update.
[1477]	93	!
[1647]	94	! VI. Surface and sub-surface soil temperature.
[1477]	95	!
[1647]	96	! VII. Perform diagnostics and write output files.
[1477]	97	!
	98	!
[253]	99	! arguments
	100	! ---------
	101	!
[1477]	102	! INPUT
[253]	103	! -----
[1477]	104	!
[253]	105	! ngrid Size of the horizontal grid.
	106	! nlayer Number of vertical layers.
[1477]	107	! nq Number of advected fields.
	108	! nametrac Name of corresponding advected fields.
	109	!
	110	! firstcall True at the first call.
	111	! lastcall True at the last call.
	112	!
	113	! pday Number of days counted from the North. Spring equinoxe.
	114	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
	115	! ptimestep timestep (s).
	116	!
	117	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
	118	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
	119	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
	120	!
	121	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
	122	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
	123	!
	124	! pt(ngrid,nlayer) Temperature (K).
	125	!
	126	! pq(ngrid,nlayer,nq) Advected fields.
	127	!
[1216]	128	! pudyn(ngrid,nlayer) \
[253]	129	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
[1477]	130	! ptdyn(ngrid,nlayer) / corresponding variables.
[253]	131	! pqdyn(ngrid,nlayer,nq) /
[1312]	132	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
[253]	133	!
[1477]	134	! OUTPUT
[253]	135	! ------
	136	!
[1308]	137	! pdu(ngrid,nlayer) \
	138	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
	139	! pdt(ngrid,nlayer) / variables due to physical processes.
	140	! pdq(ngrid,nlayer) /
[253]	141	! pdpsrf(ngrid) /
	142	!
	143	!
	144	! Authors
	145	! -------
[1524]	146	! Frederic Hourdin 15/10/93
	147	! Francois Forget 1994
	148	! Christophe Hourdin 02/1997
[253]	149	! Subroutine completely rewritten by F. Forget (01/2000)
	150	! Water ice clouds: Franck Montmessin (update 06/2003)
	151	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	152	! New correlated-k radiative scheme: R. Wordsworth (2009)
	153	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
	154	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
	155	! To F90: R. Wordsworth (2010)
[594]	156	! New turbulent diffusion scheme: J. Leconte (2012)
[716]	157	! Loops converted to F90 matrix format: J. Leconte (2012)
[787]	158	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
[1477]	159	! Purge of the code : M. Turbet (2015)
[1672]	160	! Fork for Titan : J. Vatant d'Ollone (2017)
	161	! + clean of all too-generic (ocean, water, co2 ...) routines
	162	! + Titan's chemistry
[1897]	163	! Microphysical moment model - J.Burgalat / J.Vatant d'Ollone (2017-2018)
[1647]	164	!============================================================================================
[253]	165
	166
	167	! 0. Declarations :
	168	! ------------------
	169
[1670]	170	include "netcdf.inc"
[253]	171
	172	! Arguments :
	173	! -----------
	174
[1477]	175	! INPUTS:
[253]	176	! -------
	177
[1477]	178	integer,intent(in) :: ngrid ! Number of atmospheric columns.
	179	integer,intent(in) :: nlayer ! Number of atmospheric layers.
	180	integer,intent(in) :: nq ! Number of tracers.
[1984]	181	character*30,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
[1477]	182
	183	logical,intent(in) :: firstcall ! Signals first call to physics.
	184	logical,intent(in) :: lastcall ! Signals last call to physics.
	185
	186	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
	187	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
	188	real,intent(in) :: ptimestep ! Physics timestep (s).
	189	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	190	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	191	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
	192	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
	193	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
	194	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
	195	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	196	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
[253]	197
[1477]	198	! OUTPUTS:
[253]	199	! --------
	200
[1477]	201	! Physical tendencies :
	202
	203	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
	204	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
	205	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
	206	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
	207	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
	208
[253]	209	! Local saved variables:
	210	! ----------------------
[1622]	211	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
	212	integer,save :: icount ! Counter of calls to physiq during the run.
[2366]	213	!$OMP THREADPRIVATE(day_ini,icount)
[1622]	214
[253]	215	! Local variables :
	216	! -----------------
	217
[1477]	218	real zh(ngrid,nlayer) ! Potential temperature (K).
	219	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
	220
[1670]	221	integer l,ig,ierr,iq,nw,isoil
[1161]	222
[1477]	223	! FOR DIAGNOSTIC :
[253]	224
[1477]	225	real zls ! Solar longitude (radians).
	226	real zlss ! Sub solar point longitude (radians).
	227	real zday ! Date (time since Ls=0, calculated in sols).
[2097]	228	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers (ref : local surf).
	229	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries (ref : local surf).
	230	real zzlay_eff(ngrid,nlayer) ! Effective altitude at the middle of the atmospheric layers (ref : geoid ).
	231	real zzlev_eff(ngrid,nlayer+1) ! Effective altitude at the atmospheric layer boundaries ( ref : geoid ).
	232
[1477]	233	! TENDENCIES due to various processes :
[253]	234
[1477]	235	! For Surface Temperature : (K/s)
	236	real zdtsurf(ngrid) ! Cumulated tendencies.
	237	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
	238	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
	239
	240	! For Atmospheric Temperatures : (K/s)
[1647]	241	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
[1477]	242	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
	243	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
	244
	245	! For Surface Tracers : (kg/m2/s)
	246	real dqsurf(ngrid,nq) ! Cumulated tendencies.
	247	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
	248	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
	249
	250	! For Tracers : (kg/kg_of_air/s)
	251	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
	252	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
	253	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
	254	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
[1672]	255
	256	real zdqchi(ngrid,nlayer,nq) ! Chemical tendency ( chemistry routine ).
[1795]	257
	258	real zdqmufi(ngrid,nlayer,nq) ! Microphysical tendency.
[1926]	259
	260	real zdqfibar(ngrid,nlayer,nq) ! For 2D chemistry
	261	real zdqmufibar(ngrid,nlayer,nq) ! For 2D chemistry
[1477]	262
	263	! For Winds : (m/s/s)
	264	real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
	265	real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
	266	real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	267	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	268	real zdhadj(ngrid,nlayer) ! Convadj routine.
	269
	270	! For Pressure and Mass :
	271	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
	272	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
	273	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
[253]	274
[1477]	275
	276
	277	! Local variables for LOCAL CALCULATIONS:
	278	! ---------------------------------------
[787]	279	real zflubid(ngrid)
[1308]	280	real zplanck(ngrid),zpopsk(ngrid,nlayer)
[253]	281	real ztim1,ztim2,ztim3, z1,z2
	282	real ztime_fin
[1308]	283	real zdh(ngrid,nlayer)
[1194]	284	real gmplanet
[1297]	285	real taux(ngrid),tauy(ngrid)
[1194]	286
[253]	287
[1477]	288	! local variables for DIAGNOSTICS : (diagfi & stat)
	289	! -------------------------------------------------
	290	real ps(ngrid) ! Surface Pressure.
	291	real zt(ngrid,nlayer) ! Atmospheric Temperature.
	292	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
	293	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
	294	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
	295	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
[1637]	296	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
[253]	297
[2002]	298	real zhorizwind(ngrid,nlayer) ! Horizontal Wind ( sqrt(u*+vv))
	299
[1477]	300	real vmr(ngrid,nlayer) ! volume mixing ratio
[253]	301	real time_phys
[597]	302
[1477]	303	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
	304
[594]	305	! to test energy conservation (RW+JL)
[1308]	306	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
[651]	307	real dEtot, dEtots, AtmToSurf_TurbFlux
[959]	308	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
[1315]	309	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
[1308]	310	real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
[787]	311	real dEdiffs(ngrid),dEdiff(ngrid)
[1477]	312
[594]	313	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
[1477]	314
[1295]	315	real dItot, dItot_tmp, dVtot, dVtot_tmp
[1647]	316
[1295]	317	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
[1477]	318
	319
[1482]	320	! For Clear Sky Case.
	321	real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
	322	real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux.
	323	real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation.
[1647]	324	real tf, ntf
[253]	325
[1822]	326	! Miscellaneous :
	327	character(len=10) :: tmp1
	328	character(len=10) :: tmp2
[2133]	329
	330	character*2 :: str2
[253]	331
[2366]	332	#ifndef MESOSCALE
	333
[1903]	334	! Local variables for Titan chemistry and microphysics
	335	! ----------------------------------------------------
	336
[1914]	337	real,save :: ctimestep ! Chemistry timestep (s)
	338	!$OMP THREADPRIVATE(ctimestep)
[1672]	339
[1903]	340	! Chemical tracers in molar fraction
[1914]	341	real, dimension(ngrid,nlayer,nkim) :: ychim ! (mol/mol)
[1926]	342	real, dimension(ngrid,nlayer,nkim) :: ychimbar ! For 2D chemistry
[1672]	343
[1958]	344	! Molar fraction tendencies ( chemistry, condensation and evaporation ) for tracers (mol/mol/s)
[1960]	345	real, dimension(ngrid,nlayer,nq) :: dyccond ! Condensation rate. NB : for all tracers, as we want to use indx on it.
[1958]	346	real, dimension(ngrid,nlayer,nq) :: dyccondbar ! For 2D chemistry
[1960]	347	real, dimension(ngrid) :: dycevapCH4 ! Surface "pseudo-evaporation" rate (forcing constant surface humidity).
[1672]	348
[1903]	349	! Saturation profiles
[1966]	350	real, dimension(ngrid,nlayer,nkim) :: ysat ! (mol/mol)
[1903]	351
[1789]	352
[2050]	353	real :: i2e(ngrid,nlayer) ! int 2 ext factor ( X.kg-1 -> X.m-3 for diags )
[1958]	354
[2366]	355	#ifdef USE_QTEST
[1897]	356	real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
	357	!$OMP THREADPRIVATE(tpq)
[1915]	358	real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
[2366]	359	#endif
[1897]	360
[2242]	361	logical file_ok
[1897]	362
[1672]	363	!-----------------------------------------------------------------------------
[1795]	364	! Interface to calmufi
	365	! --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module
[1897]	366	! (to have an explicit interface generated by the compiler).
[1795]	367	! Or one can put calmufi in MMP_GCM module (in muphytitan).
	368	INTERFACE
[1947]	369	SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq)
[1902]	370	REAL(kind=8), INTENT(IN) :: dt !! Physics timestep (s).
[1795]	371	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev !! Pressure levels (Pa).
	372	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev !! Altitude levels (m).
	373	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play !! Pressure layers (Pa).
	374	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay !! Altitude at the center of each layer (m).
[1947]	375	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d !! Latitude-Altitude depending gravitational acceleration (m.s-2).
[1795]	376	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp !! Temperature at the center of each layer (K).
[2242]	377	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers ($X.kg^{-1}}$).
	378	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: zdqfi !! Tendency from former processes for tracers ($X.kg^{-1}}$).
	379	REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Microphysical tendency for tracers ($X.kg^{-1}}$).
[1795]	380	END SUBROUTINE calmufi
	381	END INTERFACE
[2366]	382
	383	#endif
[1672]	384
[1477]	385	!==================================================================================================
[253]	386
	387	! -----------------
[1477]	388	! I. INITIALISATION
	389	! -----------------
[253]	390
[1477]	391	! --------------------------------
	392	! I.1 First Call Initialisation.
	393	! --------------------------------
[253]	394	if (firstcall) then
[2366]	395
	396	#ifdef USE_QTEST
[1897]	397	allocate(tpq(ngrid,nlayer,nq))
	398	tpq(:,:,:) = pq(:,:,:)
[2366]	399	#endif
[1844]	400	! Initialisation of nmicro as well as tracers names, indexes ...
[1903]	401	if (ngrid.ne.1) then ! Already done in rcm1d
[1926]	402	call initracer2(nq,nametrac) ! WARNING JB (27/03/2018): should be wrapped in an OMP SINGLE statement (see module notes)
[1844]	403	endif
[1795]	404
[2366]	405	! Allocate saved arrays (except for 1D model, where this has already been done)
	406	#ifndef MESOSCALE
	407	if (ngrid>1) call phys_state_var_init(nq)
[2291]	408	#endif
	409
[1477]	410	! Variables set to 0
[253]	411	! ~~~~~~~~~~~~~~~~~~
[2242]	412	dtrad(:,:) = 0.D0
	413	fluxrad(:) = 0.D0
	414	zdtsw(:,:) = 0.D0
	415	zdtlw(:,:) = 0.D0
[726]	416
[1822]	417	! Initialize setup for correlated-k radiative transfer
[2242]	418	! JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics.
	419	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1822]	420
	421	if (corrk) then
	422
	423	call system('rm -f surf_vals_long.out')
	424
	425	write( tmp1, '(i3)' ) L_NSPECTI
	426	write( tmp2, '(i3)' ) L_NSPECTV
	427	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	428	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	429
	430	call setspi ! Basic infrared properties.
	431	call setspv ! Basic visible properties.
	432	call sugas_corrk ! Set up gaseous absorption properties.
	433
[2133]	434	OLR_nu(:,:) = 0.D0
	435	OSR_nu(:,:) = 0.D0
[1822]	436
[2133]	437	int_dtaui(:,:,:) = 0.D0
	438	int_dtauv(:,:,:) = 0.D0
[2372]	439
	440	#ifndef MESOSCALE
[2242]	441	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
	442	haze_opt_file=trim(datadir)//'/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
	443	inquire(file=trim(haze_opt_file),exist=file_ok)
	444	if(.not.file_ok) then
	445	write(,) 'The file ',TRIM(haze_opt_file),' with the haze optical properties'
	446	write(,) 'was not found by optci.F90 ! Check in ', TRIM(datadir)
	447	write(,) 'that you have the one corresponding to the given spectral resolution !!'
	448	write(,) 'Meanwhile I abort ...'
	449	call abort
	450	endif
	451	ENDIF
[2372]	452	#endif
	453
[1822]	454	endif
	455
[2366]	456	#ifndef MESOSCALE
[1966]	457	! Initialize names and timestep for chemistry
	458	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1795]	459
[1903]	460	if ( callchim ) then
[1795]	461
[1914]	462	if ( moyzon_ch .and. ngrid.eq.1 ) then
	463	print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
	464	print *, "Please desactivate zonal mean for 1D !"
	465	stop
	466	endif
[2366]	467
[1903]	468	! Chemistry timestep
[1795]	469	ctimestep = ptimestep*REAL(ichim)
	470
[1477]	471	endif
[253]	472
[1795]	473	! Initialize microphysics.
	474	! ~~~~~~~~~~~~~~~~~~~~~~~~
[726]	475
[1795]	476	IF ( callmufi ) THEN
[2242]	477	! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
[1926]	478	call inimufi(ptimestep)
[1795]	479
[1926]	480	! initialize microphysics diagnostics arrays.
	481	call ini_diag_arrays(ngrid,nlayer,nice)
	482
[1795]	483	ENDIF
[2366]	484	#endif
[1795]	485
[1896]	486	#ifdef CPP_XIOS
	487	! Initialize XIOS context
	488	write(,) "physiq: call wxios_context_init"
	489	CALL wxios_context_init
	490	#endif
	491
[1477]	492	! Read 'startfi.nc' file.
[253]	493	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2291]	494	#ifndef MESOSCALE
[2366]	495	call phyetat0(startphy_file, &
	496	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
[1789]	497	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,tankCH4)
[2291]	498	#else
	499	emis(:)=0.0
	500	q2(:,:)=0.0
	501	qsurf(:,:)=0.0
[2366]	502	tankCH4(:)=0.0
[2291]	503	day_ini = pday
	504	#endif
	505
	506	#ifndef MESOSCALE
[1670]	507	if (.not.startphy_file) then
	508	! additionnal "academic" initialization of physics
	509	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
	510	tsurf(:)=pt(:,1)
	511	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	512	do isoil=1,nsoilmx
[1914]	513	tsoil(:,isoil)=tsurf(:)
[1670]	514	enddo
	515	if (is_master) write(,) "Physiq: initializing day_ini to pdat !"
	516	day_ini=pday
	517	endif
[2291]	518	#endif
[253]	519
	520	if (pday.ne.day_ini) then
	521	write(,) "ERROR in physiq.F90:"
	522	write(,) "bad synchronization between physics and dynamics"
	523	write(,) "dynamics day: ",pday
	524	write(,) "physics day: ",day_ini
	525	stop
	526	endif
	527	write (,) 'In physiq day_ini =', day_ini
	528
[1482]	529	! Initialize albedo calculation.
	530	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2242]	531	albedo(:,:)=0.D0
	532	albedo_bareground(:)=0.D0
[1647]	533	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
[1482]	534
	535	! Initialize orbital calculation.
	536	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	537	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	538
[2366]	539
[253]	540	if(tlocked)then
	541	print*,'Planet is tidally locked at resonance n=',nres
	542	print*,'Make sure you have the right rotation rate!!!'
	543	endif
	544
[1297]	545
[1477]	546	! Initialize soil.
	547	! ~~~~~~~~~~~~~~~~
[253]	548	if (callsoil) then
[1477]	549
[787]	550	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
[1477]	551	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[1297]	552
[1477]	553	else ! else of 'callsoil'.
	554
[253]	555	print*,'WARNING! Thermal conduction in the soil turned off'
[918]	556	capcal(:)=1.e6
[952]	557	fluxgrd(:)=intheat
	558	print*,'Flux from ground = ',intheat,' W m^-2'
[1477]	559
	560	endif ! end of 'callsoil'.
	561
[253]	562	icount=1
[1477]	563
[253]	564
[1477]	565	! Initialize surface history variable.
[253]	566	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	567	qsurf_hist(:,:)=qsurf(:,:)
[253]	568
[1637]	569	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
	570	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	571	ztprevious(:,:)=pt(:,:)
[1637]	572	zuprevious(:,:)=pu(:,:)
[253]	573
	574
[1477]	575	if(meanOLR)then
	576	call system('rm -f rad_bal.out') ! to record global radiative balance.
	577	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
	578	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
[253]	579	endif
	580
[2291]	581	#ifndef MESOSCALE
[2366]	582	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
[1542]	583	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
[2366]	584	ptimestep,pday+nday,time_phys,cell_area, &
	585	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[1904]	586	endif
[2366]	587	#endif
[1843]	588
[1622]	589	! XIOS outputs
	590	#ifdef CPP_XIOS
	591
	592	write(,) "physiq: call initialize_xios_output"
	593	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
	594	presnivs,pseudoalt)
	595	#endif
[1896]	596	write(,) "physiq: end of firstcall"
[1477]	597	endif ! end of 'firstcall'
[253]	598
[1477]	599	! ------------------------------------------------------
	600	! I.2 Initializations done at every physical timestep:
	601	! ------------------------------------------------------
	602
[1622]	603	#ifdef CPP_XIOS
	604	! update XIOS time/calendar
	605	call update_xios_timestep
	606	#endif
	607
[1477]	608	! Initialize various variables
	609	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1647]	610
[2242]	611	pdt(:,:) = 0.D0
	612	zdtsurf(:) = 0.D0
	613	pdq(:,:,:) = 0.D0
	614	dqsurf(:,:) = 0.D0
	615	pdu(:,:) = 0.D0
	616	pdv(:,:) = 0.D0
	617	pdpsrf(:) = 0.D0
	618	zflubid(:) = 0.D0
	619	taux(:) = 0.D0
	620	tauy(:) = 0.D0
[253]	621
[1477]	622	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
[1297]	623
[1477]	624	! Compute Stellar Longitude (Ls), and orbital parameters.
	625	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	626	if (season) then
	627	call stellarlong(zday,zls)
	628	else
[2366]	629	call stellarlong(noseason_day,zls)
[253]	630	end if
	631
[1329]	632	call orbite(zls,dist_star,declin,right_ascen)
[1477]	633
[1329]	634	if (tlocked) then
	635	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
	636	elseif (diurnal) then
[1524]	637	zlss=-2.pi(zday-.5)
[1329]	638	else if(diurnal .eqv. .false.) then
	639	zlss=9999.
	640	endif
[1194]	641
	642
[1477]	643	! Compute variations of g with latitude (oblate case).
	644	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	645	if (oblate .eqv. .false.) then
[1947]	646	gzlat(:,:) = g
[1477]	647	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
[1947]	648	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute gzlat (else set oblate=.false.)'
[1477]	649	call abort
	650	else
	651	gmplanet = MassPlanetgrav1e24
	652	do ig=1,ngrid
[1947]	653	gzlat(ig,:)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - latitude(ig))))
[1477]	654	end do
	655	endif
[1947]	656
	657	! Compute altitudes with the geopotential coming from the dynamics.
	658	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	659
[1947]	660	if (eff_gz .eqv. .false.) then
	661
	662	do l=1,nlayer
[2097]	663	zzlay(:,l) = pphi(:,l) / gzlat(:,l) ! Reference = local surface
[1947]	664	enddo
	665
	666	else ! In this case we consider variations of g with altitude
	667
	668	do l=1,nlayer
[2097]	669	zzlay(:,l) = gradrad / ( g*rad - ( pphi(:,l) + phisfi(:) )) - rad
[1947]	670	gzlat(:,l) = gradrad / ( rad + zzlay(:,l) )**2
	671	end do
	672
	673	endif ! if eff_gz
[728]	674
[1914]	675	zzlev(:,1)=0.
	676	zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
[2097]	677	! JVO 19 : This altitude is indeed dummy for the GCM and fits ptop=0
	678	! but for upper chemistry that's a pb -> we anyway redefine it just after ..
[728]	679
[253]	680	do l=2,nlayer
	681	do ig=1,ngrid
	682	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	683	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	684	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	685	enddo
[1477]	686	enddo
[253]	687
[2097]	688	! Effective altitudes ( eg needed for chemistry ) with correct g, and with reference to the geoid
	689	! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
	690	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2366]	691	#ifndef MESOSCALE
[2097]	692	if (moyzon_ch) then ! Zonal averages
[1672]	693
[2097]	694	zzlaybar(1,:)=gradrad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad ! reference = geoid
[1947]	695	zzlevbar(1,1)=zphisbar(1)/g
[1672]	696	DO l=2,nlayer
[1947]	697	z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplaybar(1,l-1)-zplevbar(1,l))
[1672]	698	z2=(zplevbar(1,l) +zplaybar(1,l))/(zplevbar(1,l) -zplaybar(1,l))
	699	zzlevbar(1,l)=(z1zzlaybar(1,l-1)+z2zzlaybar(1,l))/(z1+z2)
	700	ENDDO
	701	zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
	702
	703	DO ig=2,ngrid
[1947]	704	if (latitude(ig).ne.latitude(ig-1)) then
	705	DO l=1,nlayer
[2097]	706	zzlaybar(ig,l)=gradrad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
[1672]	707	ENDDO
	708	zzlevbar(ig,1)=zphisbar(ig)/g
	709	DO l=2,nlayer
[1947]	710	z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplaybar(ig,l-1)-zplevbar(ig,l))
[1672]	711	z2=(zplevbar(ig,l) +zplaybar(ig,l))/(zplevbar(ig,l) -zplaybar(ig,l))
	712	zzlevbar(ig,l)=(z1zzlaybar(ig,l-1)+z2zzlaybar(ig,l))/(z1+z2)
	713	ENDDO
[1947]	714	zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
[1672]	715	else
	716	zzlaybar(ig,:)=zzlaybar(ig-1,:)
	717	zzlevbar(ig,:)=zzlevbar(ig-1,:)
[1947]	718	endif
[1672]	719	ENDDO
	720
[2097]	721	else ! if not moyzon
[2366]	722	#endif
[2097]	723
[2098]	724	DO ig=1,ngrid
	725	DO l=1,nlayer
	726	zzlay_eff(ig,l)=gradrad/(g*rad-(pphi(ig,l)+phisfi(ig)))-rad ! reference = geoid
	727	ENDDO
	728	zzlev_eff(ig,1)=phisfi(ig)/g
	729	DO l=2,nlayer
[2097]	730	z1=(pplay(ig,l-1)+pplev(ig,l))/ (pplay(ig,l-1)-pplev(ig,l))
	731	z2=(pplev(ig,l) +pplay(ig,l))/(pplev(ig,l) -pplay(ig,l))
	732	zzlev_eff(ig,l)=(z1zzlay_eff(ig,l-1)+z2zzlay_eff(ig,l))/(z1+z2)
[2098]	733	ENDDO
	734	zzlev_eff(ig,nlayer+1)=zzlay_eff(ig,nlayer)+(zzlay_eff(ig,nlayer)-zzlev_eff(ig,nlayer))
	735	ENDDO
[2097]	736
[2366]	737	#ifndef MESOSCALE
[1672]	738	endif ! moyzon
[2366]	739	#endif
[1795]	740
[1672]	741	! -------------------------------------------------------------------------------------
[1477]	742	! Compute potential temperature
	743	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
	744	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[597]	745	do l=1,nlayer
[787]	746	do ig=1,ngrid
[253]	747	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
[597]	748	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
[1947]	749	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/gzlat(ig,l)
[1542]	750	massarea(ig,l)=mass(ig,l)*cell_area(ig)
[253]	751	enddo
	752	enddo
	753
[1312]	754	! Compute vertical velocity (m/s) from vertical mass flux
[1346]	755	! w = F / (rhoarea) and rho = P/(rT)
[1477]	756	! But first linearly interpolate mass flux to mid-layers
	757	do l=1,nlayer-1
[1914]	758	pw(:,l)=0.5*(flxw(:,l)+flxw(:,l+1))
[1477]	759	enddo
[1914]	760	pw(:,nlayer)=0.5*flxw(:,nlayer) ! since flxw(nlayer+1)=0
[1477]	761	do l=1,nlayer
[1914]	762	pw(:,l)=(pw(:,l)rpt(:,l)) / (pplay(:,l)*cell_area(:))
[1477]	763	enddo
[1194]	764
[1477]	765	!---------------------------------
	766	! II. Compute radiative tendencies
	767	!---------------------------------
[253]	768
	769	if (callrad) then
[526]	770	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[253]	771
[1477]	772	! Compute local stellar zenith angles
	773	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	774	if (tlocked) then
	775	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
	776	ztim1=SIN(declin)
	777	ztim2=COS(declin)*COS(zlss)
	778	ztim3=COS(declin)*SIN(zlss)
[253]	779
[1477]	780	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	781	ztim1,ztim2,ztim3,mu0,fract, flatten)
[253]	782
[1477]	783	elseif (diurnal) then
	784	ztim1=SIN(declin)
	785	ztim2=COS(declin)COS(2.pi*(zday-.5))
	786	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
[253]	787
[1477]	788	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	789	ztim1,ztim2,ztim3,mu0,fract, flatten)
	790	else if(diurnal .eqv. .false.) then
[253]	791
[1542]	792	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
[1161]	793	! WARNING: this function appears not to work in 1D
[253]	794
[1477]	795	endif
[1161]	796
[1477]	797	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
[1429]	798	if(rings_shadow) then
	799	call call_rings(ngrid, ptime, pday, diurnal)
	800	endif
[1133]	801
[1329]	802
[1477]	803	if (corrk) then
	804
	805	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	806	! II.a Call correlated-k radiative transfer scheme
	807	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	808
[1648]	809	call call_profilgases(nlayer)
	810
[1477]	811	! standard callcorrk
[2050]	812	call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
[1482]	813	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
[1788]	814	tsurf,fract,dist_star, &
[1482]	815	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
	816	fluxsurfabs_sw,fluxtop_lw, &
	817	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
[2133]	818	int_dtaui,int_dtauv,lastcall)
[1297]	819
[1482]	820	! Radiative flux from the sky absorbed by the surface (W.m-2).
[1477]	821	GSR=0.0
[1914]	822	fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)+fluxsurfabs_sw(:)
[253]	823
[1477]	824	!if(noradsurf)then ! no lower surface; SW flux just disappears
[1914]	825	! GSR = SUM(fluxsurf_sw(:)*cell_area(:))/totarea
	826	! fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)
[1477]	827	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
	828	!endif
[253]	829
[1477]	830	! Net atmospheric radiative heating rate (K.s-1)
[1914]	831	dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:)
[1498]	832
[1477]	833	elseif(newtonian)then
[1482]	834
	835	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	836	! II.b Call Newtonian cooling scheme
	837	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1477]	838	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
[253]	839
[1914]	840	zdtsurf(:) = +(pt(:,1)-tsurf(:))/ptimestep
[1477]	841	! e.g. surface becomes proxy for 1st atmospheric layer ?
[253]	842
[1477]	843	else
	844
	845	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	846	! II.c Atmosphere has no radiative effect
	847	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1914]	848	fluxtop_dn(:) = fract(:)mu0(:)Fat1AU/dist_star**2
[1477]	849	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	850	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	851	endif
[1914]	852	fluxsurf_sw(:) = fluxtop_dn(:)
[1482]	853	print*,'------------WARNING---WARNING------------' ! by MT2015.
	854	print*,'You are in corrk=false mode, '
[1498]	855	print*,'and the surface albedo is taken equal to the first visible spectral value'
	856
[1914]	857	fluxsurfabs_sw(:) = fluxtop_dn(:)*(1.-albedo(:,1))
	858	fluxrad_sky(:) = fluxsurfabs_sw(:)
	859	fluxtop_lw(:) = emis(:)sigmatsurf(:)**4
[253]	860
[2242]	861	dtrad(:,:)=0.D0 ! no atmospheric radiative heating
[253]	862
[1477]	863	endif ! end of corrk
[253]	864
[1477]	865	endif ! of if(mod(icount-1,iradia).eq.0)
[787]	866
[253]	867
[1477]	868	! Transformation of the radiative tendencies
	869	! ------------------------------------------
[1914]	870	zplanck(:)=tsurf(:)*tsurf(:)
	871	zplanck(:)=emis(:)sigmazplanck(:)*zplanck(:)
	872	fluxrad(:)=fluxrad_sky(:)-zplanck(:)
	873	pdt(:,:)=pdt(:,:)+dtrad(:,:)
[1477]	874
	875	! Test of energy conservation
	876	!----------------------------
[253]	877	if(enertest)then
[1524]	878	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
	879	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
[1542]	880	!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
	881	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
	882	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
[1524]	883	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	884	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
	885	if (is_master) then
[1477]	886	print*,'---------------------------------------------------------------'
	887	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	888	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	889	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	890	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	891	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	892	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
[1524]	893	endif
[1477]	894	endif ! end of 'enertest'
[253]	895
	896	endif ! of if (callrad)
	897
	898
[1477]	899
	900	! --------------------------------------------
	901	! III. Vertical diffusion (turbulent mixing) :
	902	! --------------------------------------------
	903
[253]	904	if (calldifv) then
[526]	905
[1914]	906	zflubid(:)=fluxrad(:)+fluxgrd(:)
[253]	907
[1477]	908	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
[1524]	909	if (UseTurbDiff) then
	910
[1647]	911	call turbdiff(ngrid,nlayer,nq, &
[1477]	912	ptimestep,capcal,lwrite, &
	913	pplay,pplev,zzlay,zzlev,z0, &
	914	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	915	pdt,pdq,zflubid, &
	916	zdudif,zdvdif,zdtdif,zdtsdif, &
[1647]	917	sensibFlux,q2,zdqdif,zdqsdif, &
[1477]	918	taux,tauy,lastcall)
[594]	919
[1524]	920	else
	921
[1914]	922	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
[594]	923
[1647]	924	call vdifc(ngrid,nlayer,nq,zpopsk, &
[1477]	925	ptimestep,capcal,lwrite, &
	926	pplay,pplev,zzlay,zzlev,z0, &
	927	pu,pv,zh,pq,tsurf,emis,qsurf, &
	928	zdh,pdq,zflubid, &
	929	zdudif,zdvdif,zdhdif,zdtsdif, &
[1524]	930	sensibFlux,q2,zdqdif,zdqsdif, &
[1477]	931	taux,tauy,lastcall)
[253]	932
[1914]	933	zdtdif(:,:)=zdhdif(:,:)*zpopsk(:,:) ! for diagnostic only
	934	zdqevap(:,:)=0.
[594]	935
[1477]	936	end if !end of 'UseTurbDiff'
[594]	937
[2291]	938	if (.not. turb_resolved) then
	939	pdv(:,:)=pdv(:,:)+zdvdif(:,:)
	940	pdu(:,:)=pdu(:,:)+zdudif(:,:)
	941	pdt(:,:)=pdt(:,:)+zdtdif(:,:)
	942	endif
	943	zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
[1477]	944
[253]	945	if (tracer) then
[1914]	946	pdq(:,:,:)=pdq(:,:,:)+ zdqdif(:,:,:)
	947	dqsurf(:,:)=dqsurf(:,:) + zdqsdif(:,:)
[253]	948	end if ! of if (tracer)
	949
[1477]	950
	951	! test energy conservation
[253]	952	!-------------------------
	953	if(enertest)then
[1477]	954
[1524]	955	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
[253]	956	do ig = 1, ngrid
[1524]	957	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
	958	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
[253]	959	enddo
[1477]	960
[1542]	961	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
[1524]	962	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
[1542]	963	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
	964	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
[1477]	965
[1524]	966	if (is_master) then
[1477]	967
	968	if (UseTurbDiff) then
[1524]	969	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	970	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
[1477]	971	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
[1524]	972	else
	973	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	974	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	975	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	976	end if
	977	endif ! end of 'is_master'
[1477]	978
	979	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
	980	endif ! end of 'enertest'
[253]	981
[1477]	982	else ! calldifv
[253]	983
	984	if(.not.newtonian)then
	985
[1914]	986	zdtsurf(:) = zdtsurf(:) + (fluxrad(:) + fluxgrd(:))/capcal(:)
[253]	987
	988	endif
	989
[1477]	990	endif ! end of 'calldifv'
[253]	991
	992
[1477]	993	!----------------------------------
	994	! IV. Dry convective adjustment :
	995	!----------------------------------
[253]	996
	997	if(calladj) then
	998
[1914]	999	zdh(:,:) = pdt(:,:)/zpopsk(:,:)
[2242]	1000	zduadj(:,:)=0.D0
	1001	zdvadj(:,:)=0.D0
	1002	zdhadj(:,:)=0.D0
[253]	1003
	1004
[1477]	1005	call convadj(ngrid,nlayer,nq,ptimestep, &
	1006	pplay,pplev,zpopsk, &
	1007	pu,pv,zh,pq, &
	1008	pdu,pdv,zdh,pdq, &
	1009	zduadj,zdvadj,zdhadj, &
	1010	zdqadj)
[253]	1011
[1914]	1012	pdu(:,:) = pdu(:,:) + zduadj(:,:)
	1013	pdv(:,:) = pdv(:,:) + zdvadj(:,:)
	1014	pdt(:,:) = pdt(:,:) + zdhadj(:,:)*zpopsk(:,:)
	1015	zdtadj(:,:) = zdhadj(:,:)*zpopsk(:,:) ! for diagnostic only
[1283]	1016
[253]	1017	if(tracer) then
[1914]	1018	pdq(:,:,:) = pdq(:,:,:) + zdqadj(:,:,:)
[253]	1019	end if
	1020
[1477]	1021	! Test energy conservation
[253]	1022	if(enertest)then
[1524]	1023	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
[1295]	1024	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
[253]	1025	endif
	1026
[787]	1027
[1477]	1028	endif ! end of 'calladj'
	1029
[253]	1030
[1477]	1031	!---------------------------------------------
[1647]	1032	! V. Specific parameterizations for tracers
[1477]	1033	!---------------------------------------------
[253]	1034
[1477]	1035	if (tracer) then
[1914]	1036
[2366]	1037
	1038	#ifndef MESOSCALE
	1039	!! JVO 20 : For now, no chemistry or microphysics in MESOSCALE, but why not in the future ?
	1040
[1926]	1041	! -------------------
	1042	! V.1 Microphysics
	1043	! -------------------
[1795]	1044
[1926]	1045	! JVO 05/18 : We must call microphysics before chemistry, for condensation !
	1046
	1047	if (callmufi) then
[2242]	1048
	1049	zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner)
	1050
[1926]	1051	#ifdef USE_QTEST
[2242]	1052	dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi
[1947]	1053	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
[1926]	1054	tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)ptimestep ! only manipulation of tpq->ptimestep here
	1055	#else
[2242]	1056
[2109]	1057	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
[2242]	1058
[1926]	1059	pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
[2242]	1060
	1061	! Sanity check ( way safer to be done here rather than within YAMMS )
	1062	! Important : the sanity check intentionally include the former processes tendency !
	1063	! NB : Despite this sanity check there might be still some unphysical values going through :
	1064	! - Negatives, but harmless as it will be only for the output files
	1065	! just remove them in post-proc.
	1066	! - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
	1067	! you want to compute optics from radii.
	1068	WHERE ( (pq(:,:,1)+pdq(:,:,1)ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)ptimestep < 0.D0) )
	1069	pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
	1070	pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
	1071	ENDWHERE
	1072	WHERE ( (pq(:,:,3)+pdq(:,:,3)ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)ptimestep < 0.D0) )
	1073	pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
	1074	pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
	1075	ENDWHERE
[1926]	1076	#endif
	1077
	1078	! Microphysics condensation for 2D fields to sent non-saturated fields to photochem
	1079	if ( callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0 ) then
[2117]	1080	zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
[1926]	1081	call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
[1947]	1082	gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
[2242]	1083	! TODO : Add a sanity check here !
[1926]	1084	endif
	1085
	1086	endif
	1087
	1088	! -----------------
	1089	! V.2. Chemistry
	1090	! -----------------
	1091	! NB : Must be call last ( brings fields back to an equilibrium )
	1092
	1093	! Known bug ? ( JVO 18 ) : If you try to use chimi_indx instead of iq+nmicro
	1094	! it leads to weird results / crash on dev mod ( ok in debug ) ... Why ? Idk ...
	1095
[1672]	1096	if (callchim) then
	1097
[1914]	1098	! o. Convert updated tracers to molar fraction
	1099	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1100	do iq = 1,nkim
	1101	ychim(:,:,iq) = ( pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro) ) / rat_mmol(iq+nmicro)
	1102	enddo
[1672]	1103
[1926]	1104	! JVO 05/18 : We update zonally averaged fields with condensation
	1105	! as it is compulsory to have correct photochem production. But for other
	1106	! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal
	1107	! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
	1108	if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 ) then
	1109	do iq = 1,nkim
	1110	ychimbar(:,:,iq) = zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
	1111	if ( callclouds ) then
	1112	ychimbar(:,:,iq) = ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro) )
	1113	endif
	1114	enddo
	1115	endif
	1116
[1958]	1117	! i. Condensation of the 3D tracers after the transport
	1118	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1966]	1119
[1967]	1120	call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point
[1966]	1121
[2117]	1122	dyccond(:,:,:) = 0.D0 ! Default value -> no condensation
[1966]	1123
[1903]	1124	do iq=1,nkim
[1966]	1125	where ( ychim(:,:,iq).gt.ysat(:,:,iq) ) &
	1126	dyccond(:,:,iq+nmicro) = ( -ychim(:,:,iq)+ysat(:,:,iq) ) / ptimestep
[1672]	1127	enddo
	1128
[2117]	1129	if ( callclouds ) dyccond(:,:,ices_indx) = 0.D0 ! Condensation have been calculated in the cloud microphysics
[1795]	1130
[1914]	1131	do iq=1,nkim
	1132	ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro) ! update molar ychim for following calchim
[1958]	1133
	1134	pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
[1914]	1135	enddo
	1136
	1137
[1958]	1138	! ii. 2D zonally averaged fields need to condense and evap before photochem
	1139	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1926]	1140	if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 ) then
[1958]	1141
[1967]	1142	call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
[1966]	1143
[2117]	1144	dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
[1966]	1145
[1926]	1146	do iq = 1,nkim
[1966]	1147	where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) ) &
	1148	dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
[1926]	1149	enddo
	1150
[2117]	1151	if ( callclouds ) dyccondbar(:,:,ices_indx) = 0.D0 ! Condensation have been calculated in the cloud microphysics
[1926]	1152
	1153	do iq=1,nkim
	1154	ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)
	1155	enddo
	1156
[1958]	1157	! Pseudo-evap ( forcing constant surface humidity )
	1158	do ig=1,ngrid
[2117]	1159	if ( ychimbar(ig,1,7) .lt. botCH4 ) ychimbar(ig,1,7) = botCH4
[1958]	1160	enddo
	1161
[1926]	1162	endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 )
	1163
[1958]	1164	! iii. Photochemistry ( must be call after condensation (and evap of 2D) )
	1165	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1795]	1166	if( mod(icount-1,ichim).eq.0. ) then
[1672]	1167
[1795]	1168	print *, "We enter in the chemistry ..."
[1672]	1169
[1914]	1170	if (moyzon_ch) then ! 2D zonally averaged chemistry
	1171
[1926]	1172	! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
[2097]	1173	call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar, &
[1958]	1174	zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
[1904]	1175	else ! 3D chemistry (or 1D run)
[2097]	1176	call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi, &
	1177	pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
[1914]	1178	endif ! if moyzon
[1915]	1179
[1672]	1180	endif
	1181
[1958]	1182	! iv. Surface pseudo-evaporation
	1183	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1184	do ig=1,ngrid
[2117]	1185	if ( (ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4 ) then ! +dycchi because ychim not yet updated
	1186	dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
[1958]	1187	else
[2117]	1188	dycevapCH4(ig) = 0.D0
[1958]	1189	endif
	1190	enddo
	1191
	1192	pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro)
[1956]	1193
[1958]	1194	! v. Updates and positivity check
	1195	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2242]	1196	zdqchi(:,:,:) = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step
[2117]	1197
[1903]	1198	do iq=1,nkim
[1926]	1199	zdqchi(:,:,iq+nmicro) = dycchi(:,:,iq)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
[1672]	1200
[1926]	1201	where( (pq(:,:,iq+nmicro) + ( pdq(:,:,iq+nmicro)+zdqchi(:,:,iq+nmicro) )*ptimestep ) .LT. 0.) & ! When using zonal means we set the same tendency
[2117]	1202	zdqchi(:,:,iq+nmicro) = 1.D-30 - pdq(:,:,iq+nmicro) - pq(:,:,iq+nmicro)/ptimestep ! everywhere in longitude -> could lead to negs !
[1926]	1203	enddo
	1204
[1914]	1205	pdq(:,:,:) = pdq(:,:,:) + zdqchi(:,:,:)
[1672]	1206
[1795]	1207	endif ! end of 'callchim'
	1208
[2366]	1209	#endif
	1210
[1477]	1211	! ---------------
[1672]	1212	! V.3 Updates
[1477]	1213	! ---------------
[253]	1214
[1477]	1215	! Updating Atmospheric Mass and Tracers budgets.
[728]	1216	if(mass_redistrib) then
	1217
[1914]	1218	zdmassmr(:,:) = mass(:,:) * zdqevap(:,:)
[863]	1219
	1220	do ig = 1, ngrid
[1914]	1221	zdmassmr_col(ig)=SUM(zdmassmr(ig,:))
[863]	1222	enddo
[728]	1223
[1524]	1224	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
	1225	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
	1226	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
[728]	1227
[1524]	1228	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
[1647]	1229	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
[1524]	1230	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1231	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
	1232
[1914]	1233	pdq(:,:,:) = pdq(:,:,:) + zdqmr(:,:,:)
	1234	dqsurf(:,:) = dqsurf(:,:) + zdqsurfmr(:,:)
	1235	pdt(:,:) = pdt(:,:) + zdtmr(:,:)
	1236	pdu(:,:) = pdu(:,:) + zdumr(:,:)
	1237	pdv(:,:) = pdv(:,:) + zdvmr(:,:)
	1238	pdpsrf(:) = pdpsrf(:) + zdpsrfmr(:)
	1239	zdtsurf(:) = zdtsurf(:) + zdtsurfmr(:)
[1524]	1240
	1241	endif
[728]	1242
[1477]	1243	! -----------------------------
[1672]	1244	! V.4. Surface Tracer Update
[1477]	1245	! -----------------------------
[1297]	1246
[1914]	1247	qsurf(:,:) = qsurf(:,:) + ptimestep*dqsurf(:,:)
[253]	1248
[1477]	1249	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
	1250	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
[622]	1251	qsurf_hist(:,:) = qsurf(:,:)
[253]	1252
[1477]	1253	endif! end of if 'tracer'
[253]	1254
	1255
[1477]	1256	!------------------------------------------------
[1647]	1257	! VI. Surface and sub-surface soil temperature
[1477]	1258	!------------------------------------------------
[253]	1259
[1477]	1260
	1261	! Increment surface temperature
[1297]	1262
[1914]	1263	tsurf(:)=tsurf(:)+ptimestep*zdtsurf(:)
[1477]	1264
	1265	! Compute soil temperatures and subsurface heat flux.
[253]	1266	if (callsoil) then
[787]	1267	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
[1477]	1268	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[253]	1269	endif
	1270
[1297]	1271
[1477]	1272	! Test energy conservation
[253]	1273	if(enertest)then
[1542]	1274	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
[1524]	1275	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
[253]	1276	endif
	1277
	1278
[1477]	1279	!---------------------------------------------------
[1647]	1280	! VII. Perform diagnostics and write output files
[1477]	1281	!---------------------------------------------------
	1282
	1283	! Note : For output only: the actual model integration is performed in the dynamics.
[253]	1284
[1477]	1285	! Temperature, zonal and meridional winds.
[1914]	1286	zt(:,:) = pt(:,:) + pdt(:,:)*ptimestep
	1287	zu(:,:) = pu(:,:) + pdu(:,:)*ptimestep
	1288	zv(:,:) = pv(:,:) + pdv(:,:)*ptimestep
[253]	1289
[1477]	1290	! Diagnostic.
[1914]	1291	zdtdyn(:,:) = (pt(:,:)-ztprevious(:,:)) / ptimestep
	1292	ztprevious(:,:) = zt(:,:)
[253]	1293
[1914]	1294	zdudyn(:,:) = (pu(:,:)-zuprevious(:,:)) / ptimestep
	1295	zuprevious(:,:) = zu(:,:)
[1637]	1296
[253]	1297	if(firstcall)then
[2242]	1298	zdtdyn(:,:)=0.D0
	1299	zdudyn(:,:)=0.D0
[253]	1300	endif
	1301
[2002]	1302	! Horizotal wind
	1303	zhorizwind(:,:) = sqrt( zu(:,:)zu(:,:) + zv(:,:)zv(:,:) )
	1304
[1477]	1305	! Dynamical heating diagnostic.
[253]	1306	do ig=1,ngrid
[1637]	1307	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
[253]	1308	enddo
	1309
[1477]	1310	! Tracers.
[1914]	1311	zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
[253]	1312
[1477]	1313	! Surface pressure.
[1914]	1314	ps(:) = pplev(:,1) + pdpsrf(:)*ptimestep
[253]	1315
	1316
	1317
[1477]	1318	! Surface and soil temperature information
[1542]	1319	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
[1295]	1320	call planetwide_minval(tsurf(:),Ts2)
	1321	call planetwide_maxval(tsurf(:),Ts3)
[253]	1322	if(callsoil)then
[1542]	1323	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
[1722]	1324	if (is_master) then
	1325	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	1326	print*,Ts1,Ts2,Ts3,TsS
	1327	end if
[959]	1328	else
[1722]	1329	if (is_master) then
	1330	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
[1477]	1331	print*,Ts1,Ts2,Ts3
[1524]	1332	endif
[959]	1333	end if
[253]	1334
	1335
[1477]	1336	! Check the energy balance of the simulation during the run
[253]	1337	if(corrk)then
	1338
[1542]	1339	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
	1340	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
	1341	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
	1342	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
	1343	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
[787]	1344	do ig=1,ngrid
[253]	1345	if(fluxtop_dn(ig).lt.0.0)then
	1346	print*,'fluxtop_dn has gone crazy'
	1347	print*,'fluxtop_dn=',fluxtop_dn(ig)
	1348	print*,'temp= ',pt(ig,:)
	1349	print*,'pplay= ',pplay(ig,:)
	1350	call abort
	1351	endif
	1352	end do
	1353
[787]	1354	if(ngrid.eq.1)then
[253]	1355	DYN=0.0
	1356	endif
[1524]	1357
	1358	if (is_master) then
[1477]	1359	print*,' ISR ASR OLR GND DYN [W m^-2]'
	1360	print*, ISR,ASR,OLR,GND,DYN
[1524]	1361	endif
[253]	1362
[1295]	1363	if(enertest .and. is_master)then
[651]	1364	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	1365	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	1366	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
[253]	1367	endif
	1368
[1295]	1369	if(meanOLR .and. is_master)then
[1216]	1370	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
[253]	1371	! to record global radiative balance
[588]	1372	open(92,file="rad_bal.out",form='formatted',position='append')
[651]	1373	write(92,*) zday,ISR,ASR,OLR
[253]	1374	close(92)
[588]	1375	open(93,file="tem_bal.out",form='formatted',position='append')
[1295]	1376	if(callsoil)then
[1524]	1377	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	1378	else
	1379	write(93,*) zday,Ts1,Ts2,Ts3
	1380	endif
[253]	1381	close(93)
	1382	endif
	1383	endif
	1384
[1477]	1385	endif ! end of 'corrk'
[253]	1386
[651]	1387
[1477]	1388	! Diagnostic to test radiative-convective timescales in code.
[253]	1389	if(testradtimes)then
[588]	1390	open(38,file="tau_phys.out",form='formatted',position='append')
[253]	1391	ig=1
	1392	do l=1,nlayer
	1393	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	1394	enddo
	1395	close(38)
[726]	1396	print*,'As testradtimes enabled,'
	1397	print*,'exiting physics on first call'
[253]	1398	call abort
	1399	endif
	1400
[1477]	1401
[1722]	1402	if (is_master) print,'--> Ls =',zls180./pi
[1477]	1403
	1404
	1405	!----------------------------------------------------------------------
[253]	1406	! Writing NetCDF file "RESTARTFI" at the end of the run
[1477]	1407	!----------------------------------------------------------------------
	1408
[253]	1409	! Note: 'restartfi' is stored just before dynamics are stored
	1410	! in 'restart'. Between now and the writting of 'restart',
	1411	! there will have been the itau=itau+1 instruction and
	1412	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1413	! thus we store for time=time+dtvr
	1414
	1415
	1416
[1477]	1417	if(lastcall) then
	1418	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
[305]	1419
[2291]	1420	#ifndef MESOSCALE
[1477]	1421	if (ngrid.ne.1) then
	1422	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1423
	1424	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
	1425	ptimestep,ztime_fin, &
[1789]	1426	tsurf,tsoil,emis,q2,qsurf_hist,tankCH4)
[1477]	1427	endif
[2291]	1428	#endif
	1429
[1896]	1430	endif ! end of 'lastcall'
[253]	1431
[2291]	1432	#ifndef MESOSCALE
[1477]	1433	!-----------------------------------------------------------------------------------------------------
	1434	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
	1435	!
	1436	! Note 1 : output with period "ecritphy", set in "run.def"
	1437	!
	1438	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
	1439	! but its preferable to keep all the calls in one place ...
	1440	!-----------------------------------------------------------------------------------------------------
[253]	1441
	1442
[1477]	1443	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
	1444	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
	1445	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
	1446	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
	1447	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1448	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
	1449	call writediagfi(ngrid,"temp","temperature","K",3,zt)
	1450	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
	1451	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1452	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1453	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1454	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
	1455
[965]	1456	! Subsurface temperatures
[969]	1457	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1458	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
[965]	1459
[1477]	1460
[1297]	1461
[1477]	1462	! Total energy balance diagnostics
	1463	if(callrad.and.(.not.newtonian))then
	1464
[1482]	1465	!call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
[1477]	1466	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1467	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1468	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1469
[1016]	1470	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
	1471	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
	1472	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
	1473	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
	1474	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
	1475	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
[1477]	1476
[1647]	1477
	1478	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
[1477]	1479
	1480	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
	1481
	1482	endif ! end of 'callrad'
[1524]	1483
[1477]	1484	if(enertest) then
	1485
[1524]	1486	if (calldifv) then
[1477]	1487
[1524]	1488	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
[1477]	1489	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
	1490
[1524]	1491	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	1492	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	1493	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
	1494
	1495	endif
[1477]	1496
[1524]	1497	if (corrk) then
	1498	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	1499	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
	1500	endif
[1477]	1501
	1502	endif ! end of 'enertest'
[253]	1503
[2133]	1504	! Diagnostics of optical thickness
[2138]	1505	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
[2133]	1506	if (diagdtau) then
	1507	do nw=1,L_NSPECTV
	1508	write(str2,'(i2.2)') nw
[2138]	1509	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
[2133]	1510	enddo
	1511	do nw=1,L_NSPECTI
	1512	write(str2,'(i2.2)') nw
[2138]	1513	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
[2133]	1514	enddo
	1515	endif
	1516
[1663]	1517	! Temporary inclusions for winds diagnostics.
	1518	call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
[1897]	1519	call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
[1477]	1520
	1521	! Temporary inclusions for heating diagnostics.
[1663]	1522	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	1523	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	1524	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
	1525	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
[1477]	1526
	1527	! For Debugging.
[1663]	1528	call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
[1477]	1529
[253]	1530
[1477]	1531	! Output tracers.
	1532	if (tracer) then
[1808]	1533
[1926]	1534	if (callmufi) then
	1535	! Microphysical tracers are expressed in unit/m3.
[1947]	1536	! convert X.kg-1 --> X.m-3 (whereas for optics was -> X.m-2)
	1537	i2e(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer) /(zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
[1926]	1538
[1897]	1539	#ifdef USE_QTEST
	1540	! Microphysical tracers passed through dyn+phys(except mufi)
[1947]	1541	call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
	1542	call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
	1543	call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
	1544	call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
[1897]	1545	! Microphysical tracers passed through mufi only
[1947]	1546	call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(1))*i2e)
	1547	call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(2))*i2e)
	1548	call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(3))*i2e)
	1549	call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(4))*i2e)
[1897]	1550	#else
[1947]	1551	call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
	1552	call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
	1553	call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
	1554	call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
[1897]	1555	#endif
[1926]	1556
	1557	! Microphysical diagnostics
	1558	call writediagfi(ngrid,"mmd_aer_prec","Total aerosols precipitations",'m',2,mmd_aer_prec)
	1559	call writediagfi(ngrid,"mmd_aer_s_flux","Spherical aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_s_flux)
	1560	call writediagfi(ngrid,"mmd_aer_f_flux","Fractal aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_f_flux)
	1561	call writediagfi(ngrid,"mmd_rc_sph","Spherical mode caracteristic radius",'m',3,mmd_rc_sph)
	1562	call writediagfi(ngrid,"mmd_rc_fra","Fractal mode caracteristic radius",'m',3,mmd_rc_fra)
	1563
[1795]	1564	endif ! end of 'callmufi'
[253]	1565
[1926]	1566	! Chemical tracers
	1567	if (callchim) then
	1568	do iq=1,nkim
	1569	call writediagfi(ngrid,cnames(iq),cnames(iq),'mol/mol',3,zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro))
	1570	enddo
[1960]	1571	call writediagfi(ngrid,"evapCH4","Surface CH4 pseudo-evaporation rate",'mol/mol/s',2,dycevapCH4)
[1926]	1572	endif
[1808]	1573
[1477]	1574	endif ! end of 'tracer'
[253]	1575
[1622]	1576	! XIOS outputs
	1577	#ifdef CPP_XIOS
	1578	! Send fields to XIOS: (NB these fields must also be defined as
	1579	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[1915]	1580	CALL send_xios_field("ls",zls*180./pi)
	1581	CALL send_xios_field("lss",zlss*180./pi)
	1582	CALL send_xios_field("RA",right_ascen*180./pi)
	1583	CALL send_xios_field("Declin",declin*180./pi)
[1626]	1584
[1915]	1585	! Total energy balance diagnostics
	1586	if (callrad.and.(.not.newtonian)) then
	1587	CALL send_xios_field("ISR_TOA",fluxtop_dn)
	1588	CALL send_xios_field("OLR_TOA",fluxtop_lw)
	1589	endif
	1590
	1591	CALL send_xios_field("area",cell_area)
	1592	CALL send_xios_field("pphi",pphi)
[1957]	1593	CALL send_xios_field("pphis",phisfi)
[1915]	1594
[1622]	1595	CALL send_xios_field("ps",ps)
[1915]	1596	CALL send_xios_field("tsurf",tsurf)
[1622]	1597
[1957]	1598	if(enertest) then
	1599	if (calldifv) then
	1600	CALL send_xios_field("sensibFlux",sensibFlux)
	1601	endif
	1602	endif
	1603
[1915]	1604	CALL send_xios_field("temp",zt)
	1605	CALL send_xios_field("teta",zh)
[1622]	1606	CALL send_xios_field("u",zu)
	1607	CALL send_xios_field("v",zv)
[1915]	1608	CALL send_xios_field("w",pw)
	1609	CALL send_xios_field("p",pplay)
	1610
	1611	! Winds diagnostics.
	1612	CALL send_xios_field("dudif",zdudif)
	1613	CALL send_xios_field("dudyn",zdudyn)
[1896]	1614
[2002]	1615	CALL send_xios_field("horizwind",zhorizwind)
	1616
[1915]	1617	! Heating diagnostics.
	1618	CALL send_xios_field("dtsw",zdtsw)
	1619	CALL send_xios_field("dtlw",zdtlw)
	1620	CALL send_xios_field("dtrad",dtrad)
	1621	CALL send_xios_field("dtdyn",zdtdyn)
[2097]	1622	CALL send_xios_field("dtdif",zdtdif)
[1926]	1623
	1624	! Chemical tracers
	1625	if (callchim) then
[1943]	1626
	1627	! Advected fields
[1926]	1628	do iq=1,nkim
[1943]	1629	CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
[1926]	1630	enddo
[1943]	1631
	1632	! Upper chemistry fields
	1633	do iq=1,nkim
	1634	CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
	1635	enddo
	1636
	1637	! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
	1638	do iq=1,nkim
	1639
	1640	! GCM levels
	1641	do l=1,nlayer
	1642	ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
	1643	enddo
	1644	! Upper levels
	1645	do l=1,nlaykim_up
	1646	ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
	1647	enddo
	1648
	1649	CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
	1650
	1651	enddo
[1947]	1652
	1653	! Condensation tendencies ( mol/mol/s )
[1958]	1654	CALL send_xios_field("dqcond_CH4",dyccond(:,:,7+nmicro))
	1655	CALL send_xios_field("dqcond_C2H2",dyccond(:,:,10+nmicro))
	1656	CALL send_xios_field("dqcond_C2H4",dyccond(:,:,12+nmicro))
	1657	CALL send_xios_field("dqcond_C2H6",dyccond(:,:,14+nmicro))
	1658	CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
	1659	CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
[2242]	1660	CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
	1661	CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
	1662	CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
	1663	CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
	1664	CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
	1665	CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
	1666	CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
	1667	CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
	1668	CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
	1669	CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
	1670	CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
[1947]	1671
[1958]	1672	! Pseudo-evaporation flux (mol/mol/s)
[1960]	1673	CALL send_xios_field("evapCH4",dycevapCH4(:))
[1958]	1674
[1943]	1675	endif ! of 'if callchim'
[1926]	1676
	1677	! Microphysical tracers
	1678	if (callmufi) then
[1947]	1679	CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
	1680	CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
	1681	CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
	1682	CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
[1926]	1683	endif
	1684
[1896]	1685	if (lastcall.and.is_omp_master) then
	1686	write(,) "physiq: call xios_context_finalize"
	1687	call xios_context_finalize
	1688	endif
[1622]	1689	#endif
[2291]	1690	#else
	1691	!MESOSCALE outputs
	1692	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
	1693	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
	1694	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
	1695	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
	1696	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
	1697	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
	1698	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
	1699	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
	1700	#endif
[1622]	1701
[253]	1702	icount=icount+1
	1703
[751]	1704	end subroutine physiq
[1549]	1705
	1706	end module physiq_mod

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

Download in other formats:

Original Format