Context Navigation

physiq_mod.F90 @ 2509

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

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

Download in other formats:

Original Format