Context Navigation

physiq_mod.F90 @ 2436

Last change on this file since 2436 was 2428, checked in by emillour, 5 years ago
Generci GCM: Bug fix on mass_redistribution; argument rnat should be real, not integer. Turned it into a mass_redistribution_mod module. EM + YJ
Property svn:executable set to ``*
File size: 100.3 KB

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

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

Download in other formats:

Original Format