Context Navigation

physiq_mod.F90 @ 3184

Last change on this file since 3184 was 3184, checked in by afalco, 10 months ago
Pluto PCM: Added LMDZ.PLUTO, a copy of the generic model, cleaned from some unnecessary modules (water, ...) AF
Property svn:executable set to ``*
File size: 75.9 KB

Line
1	module physiq_mod
2
3	implicit none
4
5	contains
6
7	subroutine physiq(ngrid,nlayer,nq, &
8	firstcall,lastcall, &
9	pday,ptime,ptimestep, &
10	pplev,pplay,pphi, &
11	pu,pv,pt,pq, &
12	flxw, &
13	pdu,pdv,pdt,pdq,pdpsrf)
14
15	!!
16	use write_field_phy, only: Writefield_phy
17	!!
18	use ioipsl_getin_p_mod, only: getin_p
19	use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
20	use generic_cloud_common_h, only : epsi_generic, Psat_generic
21	use thermcell_mod, only: init_thermcell_mod
22	use gases_h, only: gnom, gfrac
23	use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
24	use suaer_corrk_mod, only: suaer_corrk
25	use radii_mod, only: n2_reffrad
26	use aerosol_mod, only: iniaerosol, iaero_n2
27	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
28	dryness
29	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
30	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
31	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
32	use geometry_mod, only: latitude, longitude, cell_area
33	USE comgeomfi_h, only: totarea, totarea_planet
34	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
35	alpha_lift, alpha_devil, qextrhor, &
36	igcm_dustbin, &
37	igcm_n2_ice, nesp, is_chim, is_condensable,constants_epsi_generic
38	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
39	use phyetat0_mod, only: phyetat0
40	use wstats_mod, only: callstats, wstats, mkstats
41	use phyredem, only: physdem0, physdem1
42	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
43	use mod_phys_lmdz_para, only : is_master
44	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
45	obliquit, nres, z0
46	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
47	use time_phylmdz_mod, only: daysec
48	use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
49	callrad, callsoil, nosurf, &
50	n2cond, corrk, diagdtau, &
51	diurnal, enertest, fat1au, &
52	icetstep, intheat, iradia, kastprof, &
53	lwrite, mass_redistrib, meanOLR, &
54	nearn2cond, noseason_day, &
55	season, sedimentation,generic_condensation, &
56	specOLR, &
57	startphy_file, testradtimes, &
58	tracer, UseTurbDiff, &
59	global1d, szangle
60	use generic_tracer_index_mod, only: generic_tracer_index
61	use check_fields_mod, only: check_physics_fields
62	use conc_mod, only: rnew, cpnew, ini_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	use condensation_generic_mod, only: condensation_generic
70	use datafile_mod, only: datadir
71	#ifndef MESOSCALE
72	use vertical_layers_mod, only: presnivs, pseudoalt
73	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
74	#else
75	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
76	comm_CLOUDFRAC,comm_TOTCLOUDFRAC, comm_RH, &
77	comm_HR_DYN, &
78	comm_DQICE,comm_DQVAP,comm_ALBEQ, &
79	comm_FLUXTOP_DN,comm_FLUXABS_SW, &
80	comm_FLUXTOP_LW,comm_FLUXSURF_SW, &
81	comm_FLUXSURF_LW,comm_FLXGRD, &
82	comm_DTRAIN,comm_DTLSC, &
83	comm_LATENT_HF
84
85	#endif
86
87	#ifdef CPP_XIOS
88	use xios_output_mod, only: initialize_xios_output, &
89	update_xios_timestep, &
90	send_xios_field
91	use wxios, only: wxios_context_init, xios_context_finalize
92	#endif
93
94	implicit none
95
96
97	!==================================================================
98	!
99	! Purpose
100	! -------
101	! Central subroutine for all the physics parameterisations in the
102	! universal model. Originally adapted from the Mars LMDZ model.
103	!
104	! The model can be run without or with tracer transport
105	! depending on the value of "tracer" in file "callphys.def".
106	!
107	!
108	! It includes:
109	!
110	! I. Initialization :
111	! I.1 Firstcall initializations.
112	! I.2 Initialization for every call to physiq.
113	!
114	! II. Compute radiative transfer tendencies (longwave and shortwave) :
115	! II.a Option 1 : Call correlated-k radiative transfer scheme.
116	! II.b Option 2 : Call Newtonian cooling scheme.
117	! II.c Option 3 : Atmosphere has no radiative effect.
118	!
119	! III. Vertical diffusion (turbulent mixing) :
120	!
121	! IV. Convection :
122	! IV.a Thermal plume model
123	! IV.b Dry convective adjusment
124	!
125	! V. Condensation and sublimation of gases (currently just N2).
126	!
127	! VI. Tracers
128	! VI.1. Water and water ice.
129	! VI.2 Photochemistry
130	! VI.3. Aerosols and particles.
131	! VI.4. Updates (pressure variations, surface budget).
132	! VI.5. Slab Ocean.
133	! VI.6. Surface Tracer Update.
134	!
135	! VII. Surface and sub-surface soil temperature.
136	!
137	! VIII. Perform diagnostics and write output files.
138	!
139	!
140	! arguments
141	! ---------
142	!
143	! INPUT
144	! -----
145	!
146	! ngrid Size of the horizontal grid.
147	! nlayer Number of vertical layers.
148	! nq Number of advected fields.
149	!
150	! firstcall True at the first call.
151	! lastcall True at the last call.
152	!
153	! pday Number of days counted from the North. Spring equinoxe.
154	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
155	! ptimestep timestep (s).
156	!
157	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
158	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
159	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
160	!
161	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
162	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
163	!
164	! pt(ngrid,nlayer) Temperature (K).
165	!
166	! pq(ngrid,nlayer,nq) Advected fields.
167	!
168	! pudyn(ngrid,nlayer) \
169	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
170	! ptdyn(ngrid,nlayer) / corresponding variables.
171	! pqdyn(ngrid,nlayer,nq) /
172	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
173	!
174	! OUTPUT
175	! ------
176	!
177	! pdu(ngrid,nlayer) \
178	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
179	! pdt(ngrid,nlayer) / variables due to physical processes.
180	! pdq(ngrid,nlayer) /
181	! pdpsrf(ngrid) /
182	!
183	!
184	! Authors
185	! -------
186	! Frederic Hourdin 15/10/93
187	! Francois Forget 1994
188	! Christophe Hourdin 02/1997
189	! Subroutine completely rewritten by F. Forget (01/2000)
190	! Water ice clouds: Franck Montmessin (update 06/2003)
191	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
192	! New correlated-k radiative scheme: R. Wordsworth (2009)
193	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
194	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
195	! To F90: R. Wordsworth (2010)
196	! New turbulent diffusion scheme: J. Leconte (2012)
197	! Loops converted to F90 matrix format: J. Leconte (2012)
198	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
199	! Purge of the code : M. Turbet (2015)
200	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
201	! Purge for Pluto model : A. Falco (2024)
202	!==================================================================
203
204
205	! 0. Declarations :
206	! ------------------
207
208	include "netcdf.inc"
209
210	! Arguments :
211	! -----------
212
213	! INPUTS:
214	! -------
215
216	integer,intent(in) :: ngrid ! Number of atmospheric columns.
217	integer,intent(in) :: nlayer ! Number of atmospheric layers.
218	integer,intent(in) :: nq ! Number of tracers.
219
220	logical,intent(in) :: firstcall ! Signals first call to physics.
221	logical,intent(in) :: lastcall ! Signals last call to physics.
222
223	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
224	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
225	real,intent(in) :: ptimestep ! Physics timestep (s).
226	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
227	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
228	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
229	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
230	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
231	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
232	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
233	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
234
235	! OUTPUTS:
236	! --------
237
238	! Physical tendencies :
239
240	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
241	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
242	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
243	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
244	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
245
246	! Local saved variables:
247	! ----------------------
248	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
249	integer,save :: icount ! Counter of calls to physiq during the run.
250	!$OMP THREADPRIVATE(day_ini,icount)
251
252	! Local variables :
253	! -----------------
254
255	! Aerosol (dust or ice) extinction optical depth at reference wavelength
256	! for the "naerkind" optically active aerosols:
257
258	real,save,allocatable :: aerosol(:,:,:) ! Aerosols
259	!$OMP THREADPRIVATE(aerosol)
260	real zh(ngrid,nlayer) ! Potential temperature (K).
261	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
262	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
263
264	integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
265	logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
266
267	real zls ! Solar longitude (radians).
268	real zlss ! Sub solar point longitude (radians).
269	real zday ! Date (time since Ls=0, calculated in sols).
270	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
271	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
272
273	! VARIABLES for the thermal plume model (AF24: deleted)
274
275	! TENDENCIES due to various processes :
276
277	! For Surface Temperature : (K/s)
278	real zdtsurf(ngrid) ! Cumulated tendencies.
279	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
280	real zdtsurfc(ngrid) ! Condense_n2 routine.
281	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
282	! real zdtsurf_hyd(ngrid) ! Hydrol routine.
283
284	! For Atmospheric Temperatures : (K/s)
285	real dtlscale(ngrid,nlayer) ! Largescale routine.
286	real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine.
287	real zdtc(ngrid,nlayer) ! Condense_n2 routine.
288	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
289	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
290	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
291	real zdtchim(ngrid,nlayer) ! Calchim routine.
292
293	! For Surface Tracers : (kg/m2/s)
294	real dqsurf(ngrid,nq) ! Cumulated tendencies.
295	real zdqsurfc(ngrid) ! Condense_n2 routine.
296	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
297	real zdqssed(ngrid,nq) ! Callsedim routine.
298	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
299
300	! For Tracers : (kg/kg_of_air/s)
301	real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine.
302	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
303	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
304	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
305	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
306	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
307	real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
308	real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
309	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
310	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
311	!$OMP THREADPRIVATE(zdqchim,zdqschim)
312
313	REAL array_zero1(ngrid)
314	REAL array_zero2(ngrid,nlayer)
315
316	! For Winds : (m/s/s)
317	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
318	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
319	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
320	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
321	real zdhadj(ngrid,nlayer) ! Convadj routine.
322
323	! For Pressure and Mass :
324	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
325	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
326	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
327
328
329
330	! Local variables for LOCAL CALCULATIONS:
331	! ---------------------------------------
332	real zflubid(ngrid)
333	real zplanck(ngrid),zpopsk(ngrid,nlayer)
334	real ztim1,ztim2,ztim3, z1,z2
335	real ztime_fin
336	real zdh(ngrid,nlayer)
337	real gmplanet
338	real taux(ngrid),tauy(ngrid)
339
340
341	! local variables for DIAGNOSTICS : (diagfi & stat)
342	! -------------------------------------------------
343	real ps(ngrid) ! Surface Pressure.
344	real zt(ngrid,nlayer) ! Atmospheric Temperature.
345	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
346	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
347	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
348	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
349	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
350
351	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
352	real vmr(ngrid,nlayer) ! volume mixing ratio
353	real time_phys
354
355	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
356
357	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
358
359	! included by RW for H2O Manabe scheme
360	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
361	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
362
363
364	! to test energy conservation (RW+JL)
365	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
366	real dEtot, dEtots, AtmToSurf_TurbFlux
367	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
368	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
369
370	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
371
372	real dtmoist_max,dtmoist_min
373	real dItot, dItot_tmp, dVtot, dVtot_tmp
374
375
376	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
377
378	real psat_tmp ! AF24: to remove?
379
380	real qsat_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) specific concentration at saturation (kg/kg_of_air).
381	real RH_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) Relative humidity.
382	real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation).
383	real psat_tmp_generic
384	real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
385	!$OMP THREADPRIVATE(metallicity)
386
387	real reffrad_generic_zeros_for_wrf(ngrid,nlayer) ! !!! this is temporary, it is only a list of zeros, it will be replaced when a generic aerosol will be implemented
388
389	! For Clear Sky Case. (AF24: deleted)
390
391	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
392
393	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
394
395	real,save,allocatable :: reffcol(:,:)
396	!$OMP THREADPRIVATE(reffcol)
397
398	! Non-oro GW tendencies
399	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
400	REAL d_t_hin(ngrid,nlayer)
401	! Diagnostics 2D of gw_nonoro
402	REAL zustrhi(ngrid), zvstrhi(ngrid)
403
404
405	real :: tsurf2(ngrid)
406	!! real :: flux_o(ngrid),flux_g(ngrid)
407	real :: flux_g(ngrid)
408	real :: flux_sens_lat(ngrid)
409	real :: qsurfint(ngrid,nq)
410	#ifdef MESOSCALE
411	REAL :: lsf_dt(nlayer)
412	REAL :: lsf_dq(nlayer)
413	#endif
414
415	! flags to trigger extra sanity checks
416	logical, save :: check_physics_inputs=.false.
417	logical, save :: check_physics_outputs=.false.
418	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
419
420	! Misc
421	character*2 :: str2
422	character(len=10) :: tmp1
423	character(len=10) :: tmp2
424	!==================================================================================================
425
426	! -----------------
427	! I. INITIALISATION
428	! -----------------
429
430	! --------------------------------
431	! I.1 First Call Initialisation.
432	! --------------------------------
433	if (firstcall) then
434	call getin_p("check_physics_inputs", check_physics_inputs)
435	call getin_p("check_physics_outputs", check_physics_outputs)
436
437	! Allocate saved arrays (except for 1D model, where this has already
438	! been done)
439	#ifndef MESOSCALE
440	if (ngrid>1) call phys_state_var_init(nq)
441	#endif
442
443	! Variables set to 0
444	! ~~~~~~~~~~~~~~~~~~
445	dtrad(:,:) = 0.0
446	fluxrad(:) = 0.0
447	tau_col(:) = 0.0
448	zdtsw(:,:) = 0.0
449	zdtlw(:,:) = 0.0
450
451	! Initialize tracer names, indexes and properties.
452	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
454	if (tracer) then
455	call initracer(ngrid,nq)
456	! if(photochem) then !AF24: removed
457	endif
458	! Initialize aerosol indexes.
459	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
460	call iniaerosol
461	! allocate related local arrays
462	! (need be allocated instead of automatic because of "naerkind")
463	allocate(aerosol(ngrid,nlayer,naerkind))
464	allocate(reffcol(ngrid,naerkind))
465
466	#ifdef CPP_XIOS
467	! Initialize XIOS context
468	write(,) "physiq: call wxios_context_init"
469	CALL wxios_context_init
470	#endif
471
472	! Read 'startfi.nc' file.
473	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
474	#ifndef MESOSCALE
475	call phyetat0(startphy_file, &
476	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
477	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
478
479	#else
480
481	day_ini = pday
482	#endif
483
484	#ifndef MESOSCALE
485	if (.not.startphy_file) then
486	! additionnal "academic" initialization of physics
487	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
488	tsurf(:)=pt(:,1)
489	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
490	do isoil=1,nsoilmx
491	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
492	enddo
493	if (is_master) write(,) "Physiq: initializing day_ini to pdat !"
494	day_ini=pday
495	endif
496	#endif
497	if (pday.ne.day_ini) then
498	write(,) "ERROR in physiq.F90:"
499	write(,) "bad synchronization between physics and dynamics"
500	write(,) "dynamics day: ",pday
501	write(,) "physics day: ",day_ini
502	stop
503	endif
504
505	write (,) 'In physiq day_ini =', day_ini
506
507	! Initialize albedo calculation.
508	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
509	albedo(:,:)=0.0
510	albedo_bareground(:)=0.0
511	albedo_snow_SPECTV(:)=0.0
512	albedo_n2_ice_SPECTV(:)=0.0
513	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
514
515	! Initialize orbital calculation.
516	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
517	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
518
519
520	! Initialize soil.
521	! ~~~~~~~~~~~~~~~~
522	if (callsoil) then
523	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
524	ptimestep,tsurf,tsoil,capcal,fluxgrd)
525	else ! else of 'callsoil'.
526	print*,'WARNING! Thermal conduction in the soil turned off'
527	capcal(:)=1.e6
528	fluxgrd(:)=intheat
529	print*,'Flux from ground = ',intheat,' W m^-2'
530	endif ! end of 'callsoil'.
531
532	icount=1
533
534	! Initialize surface history variable.
535	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
536	qsurf_hist(:,:)=qsurf(:,:)
537
538	!! call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
539
540	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
541	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
542	ztprevious(:,:)=pt(:,:)
543	zuprevious(:,:)=pu(:,:)
544
545	! Set temperature just above condensation temperature (for Early Mars)
546	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
547	if(nearn2cond) then
548	write(,)' WARNING! Starting at Tcond+1K'
549	do l=1, nlayer
550	do ig=1,ngrid
551	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
552	-pt(ig,l)) / ptimestep
553	enddo
554	enddo
555	endif
556
557	if(meanOLR)then
558	call system('rm -f rad_bal.out') ! to record global radiative balance.
559	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
560	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
561	endif
562
563
564	!! Initialize variables for water cycle ! AF24: removed
565
566	! Set metallicity for GCS
567	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
568	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
569	call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
570
571	! Set some parameters for the thermal plume model !AF24: removed
572	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
573
574	#ifndef MESOSCALE
575	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
576	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
577	ptimestep,pday+nday,time_phys,cell_area, &
578	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
579	endif
580
581	!! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
582	#endif
583	if (corrk) then
584	! We initialise the spectral grid here instead of
585	! at firstcall of callcorrk so we can output XspecIR, XspecVI
586	! when using Dynamico
587	print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
588	call getin_p("corrkdir",corrkdir)
589	print*,"corrkdir = ", corrkdir
590	write (tmp1, '(i4)') L_NSPECTI
591	write (tmp2, '(i4)') L_NSPECTV
592	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
593	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
594	call setspi !Basic infrared properties.
595	call setspv ! Basic visible properties.
596	call sugas_corrk ! Set up gaseous absorption properties.
597	call suaer_corrk ! Set up aerosol optical properties.
598	endif
599
600	!! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
601	! XIOS outputs
602	#ifdef CPP_XIOS
603
604	write(,) "physiq: call initialize_xios_output"
605	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
606	year_day,presnivs,pseudoalt,WNOI,WNOV)
607	#endif
608
609	!! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
610
611	write(,) "physiq: end of firstcall"
612	endif ! end of 'firstcall'
613
614	!! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
615	!! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
616
617	if (check_physics_inputs) then
618	!check the validity of input fields coming from the dynamics
619	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
620	endif
621
622	! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
623	! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
624
625	! ------------------------------------------------------
626	! I.2 Initializations done at every physical timestep:
627	! ------------------------------------------------------
628
629	#ifdef CPP_XIOS
630	! update XIOS time/calendar
631	call update_xios_timestep
632	#endif
633
634	! Initialize various variables
635	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
636
637	if ( .not.nearn2cond ) then
638	pdt(1:ngrid,1:nlayer) = 0.0
639	endif
640	zdtsurf(1:ngrid) = 0.0
641	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
642	dqsurf(1:ngrid,1:nq)= 0.0
643	pdu(1:ngrid,1:nlayer) = 0.0
644	pdv(1:ngrid,1:nlayer) = 0.0
645	pdpsrf(1:ngrid) = 0.0
646	zflubid(1:ngrid) = 0.0
647	flux_sens_lat(1:ngrid) = 0.0
648	taux(1:ngrid) = 0.0
649	tauy(1:ngrid) = 0.0
650
651	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
652
653	! Compute Stellar Longitude (Ls), and orbital parameters.
654	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
655	if (season) then
656	call stellarlong(zday,zls)
657	else
658	call stellarlong(noseason_day,zls)
659	end if
660
661	call orbite(zls,dist_star,declin,right_ascen)
662
663	if (diurnal) then
664	zlss=-2.pi(zday-.5)
665	else if(diurnal .eqv. .false.) then
666	zlss=9999.
667	endif
668
669
670	glat(:) = g !AF24: removed oblateness
671
672	! Compute geopotential between layers.
673	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
674	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
675	do l=1,nlayer
676	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
677	enddo
678
679	zzlev(1:ngrid,1)=0.
680
681	do l=2,nlayer
682	do ig=1,ngrid
683	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
684	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
685	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
686	enddo
687	enddo
688
689	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
690
691	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
692
693	! Compute potential temperature
694	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
695	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
696	do l=1,nlayer
697	do ig=1,ngrid
698	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
699	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
700	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
701	massarea(ig,l)=mass(ig,l)*cell_area(ig)
702	enddo
703	enddo
704
705	! Compute vertical velocity (m/s) from vertical mass flux
706	! w = F / (rhoarea) and rho = P/(rT)
707	! But first linearly interpolate mass flux to mid-layers
708	do l=1,nlayer-1
709	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
710	enddo
711	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
712	do l=1,nlayer
713	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
714	(pplay(1:ngrid,l)*cell_area(1:ngrid))
715	enddo
716	! omega in Pa/s
717	do l=1,nlayer-1
718	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
719	enddo
720	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
721	do l=1,nlayer
722	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
723	enddo
724
725	!---------------------------------
726	! II. Compute radiative tendencies
727	!---------------------------------
728
729	! Compute local stellar zenith angles
730	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
731	if (diurnal) then
732	ztim1=SIN(declin)
733	ztim2=COS(declin)COS(2.pi*(zday-.5))
734	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
735
736	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
737	ztim1,ztim2,ztim3,mu0,fract)
738	else if(diurnal .eqv. .false.) then
739
740	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
741	! WARNING: this function appears not to work in 1D
742
743	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
744	mu0 = cos(pi*szangle/180.0)
745	endif
746
747	endif
748
749	if (callrad) then
750	if( mod(icount-1,iradia).eq.0.or.lastcall) then
751
752	! Eclipse incoming sunlight !AF24: removed
753
754	!! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
755	!! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
756
757
758	if (corrk) then
759
760	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761	! II.a Call correlated-k radiative transfer scheme
762	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
763	if(kastprof)then
764	print*,'kastprof should not = true here'
765	call abort
766	endif
767	! if(water) then !AF24: removed
768
769	if(generic_condensation) then
770	do iq=1,nq
771
772	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
773
774	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
775
776	epsi_generic=constants_epsi_generic(iq)
777
778	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
779	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
780
781	endif
782	end do ! do iq=1,nq loop on tracers
783	! take into account generic condensable specie (GCS) effect on mean molecular weight
784
785	else
786	muvar(1:ngrid,1:nlayer+1)=mugaz
787	endif
788
789	! if(ok_slab_ocean) then !AF24: removed
790
791	! standard callcorrk
792	! clearsky=.false.
793	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
794	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
795	tsurf,fract,dist_star,aerosol,muvar, &
796	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
797	fluxsurfabs_sw,fluxtop_lw, &
798	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
799	int_dtaui,int_dtauv, &
800	tau_col,cloudfrac,totcloudfrac, &
801	.false.,firstcall,lastcall)
802
803	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
804	!the deep atmosphere press > 0.1 bar. Suggested by MT.
805	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
806
807	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
808	!coeff_rad=0.5
809	!coeff_rad=1.
810	!do l=1, nlayer
811	! do ig=1,ngrid
812	! if(pplay(ig,l).ge.1.d4) then
813	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
814	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
815	! endif
816	! enddo
817	!enddo
818
819	! AF24: removed CLFvarying Option
820
821	! if(ok_slab_ocean) then
822	! tsurf(:)=tsurf2(:)
823	! endif
824
825
826	! Radiative flux from the sky absorbed by the surface (W.m-2).
827	GSR=0.0
828	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
829
830	!if(noradsurf)then ! no lower surface; SW flux just disappears
831	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
832	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
833	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
834	!endif
835
836	! Net atmospheric radiative heating rate (K.s-1)
837	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
838
839	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
840	if (firstcall .and. albedo_spectral_mode) then
841	call spectral_albedo_calc(albedo_snow_SPECTV)
842	endif
843
844	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
845	! II.b Call Newtonian cooling scheme !AF24: removed
846	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
847	else
848	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
849	! II.c Atmosphere has no radiative effect
850	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
851	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
852	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
853	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
854	endif
855	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
856	print*,'------------WARNING---WARNING------------' ! by MT2015.
857	print*,'You are in corrk=false mode, '
858	print*,'and the surface albedo is taken equal to the first visible spectral value'
859
860	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
861	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
862	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
863
864	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
865
866	endif ! end of corrk
867
868	endif ! of if(mod(icount-1,iradia).eq.0)
869
870
871	! Transformation of the radiative tendencies
872	! ------------------------------------------
873	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
874	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
875	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
876	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
877
878	! Test of energy conservation
879	!----------------------------
880	if(enertest)then
881	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
882	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
883	!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
884	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
885	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
886	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
887	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
888	if (is_master) then
889	print*,'---------------------------------------------------------------'
890	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
891	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
892	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
893	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
894	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
895	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
896	endif
897	endif ! end of 'enertest'
898
899	endif ! of if (callrad)
900
901	!! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
902	!! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
903
904
905	! --------------------------------------------
906	! III. Vertical diffusion (turbulent mixing) :
907	! --------------------------------------------
908
909	if (calldifv) then
910
911	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
912
913	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
914	if (UseTurbDiff) then
915
916	call turbdiff(ngrid,nlayer,nq, &
917	ptimestep,capcal, &
918	pplay,pplev,zzlay,zzlev,z0, &
919	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
920	pdt,pdq,zflubid, &
921	zdudif,zdvdif,zdtdif,zdtsdif, &
922	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
923	taux,tauy)
924
925	else
926
927	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
928
929	call vdifc(ngrid,nlayer,nq,zpopsk, &
930	ptimestep,capcal,lwrite, &
931	pplay,pplev,zzlay,zzlev,z0, &
932	pu,pv,zh,pq,tsurf,emis,qsurf, &
933	zdh,pdq,zflubid, &
934	zdudif,zdvdif,zdhdif,zdtsdif, &
935	sensibFlux,q2,zdqdif,zdqsdif)
936
937	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
938	zdqevap(1:ngrid,1:nlayer)=0.
939
940	end if !end of 'UseTurbDiff'
941
942	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
943
944	!!! this is always done, except for turbulence-resolving simulations
945	if (.not. turb_resolved) then
946	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
947	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
948	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
949	endif
950
951	! if(ok_slab_ocean)then !AF24: removed
952	! flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
953	! endif
954
955	!! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
956
957	if (tracer) then
958	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
959	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
960	end if ! of if (tracer)
961
962	!! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
963	!! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
964
965	! test energy conservation
966	!-------------------------
967	if(enertest)then
968
969	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
970	do ig = 1, ngrid
971	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
972	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
973	enddo
974
975	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
976	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
977	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
978	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
979
980	if (is_master) then
981
982	if (UseTurbDiff) then
983	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
984	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
985	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
986	else
987	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
988	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
989	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
990	end if
991	endif ! end of 'is_master'
992
993	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
994	endif ! end of 'enertest'
995
996
997	! ! Test water conservation. !AF24: removed
998
999	else ! calldifv
1000
1001	! if(.not.newtonian)then
1002	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1003
1004	endif ! end of 'calldifv'
1005
1006
1007	!-------------------
1008	! IV. Convection :
1009	!-------------------
1010
1011	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1012	! IV.a Thermal plume model : AF24: removed
1013	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1014
1015	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1016	! IV.b Dry convective adjustment :
1017	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1018
1019	if(calladj) then
1020
1021	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1022	zduadj(1:ngrid,1:nlayer)=0.0
1023	zdvadj(1:ngrid,1:nlayer)=0.0
1024	zdhadj(1:ngrid,1:nlayer)=0.0
1025
1026
1027	call convadj(ngrid,nlayer,nq,ptimestep, &
1028	pplay,pplev,zpopsk, &
1029	pu,pv,zh,pq, &
1030	pdu,pdv,zdh,pdq, &
1031	zduadj,zdvadj,zdhadj, &
1032	zdqadj)
1033
1034	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1035	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1036	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1037	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1038
1039	if(tracer) then
1040	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1041	end if
1042
1043	! Test energy conservation
1044	if(enertest)then
1045	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1046	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1047	endif
1048
1049	! ! Test water conservation !AF24: removed
1050
1051	endif ! end of 'calladj'
1052	!----------------------------------------------
1053	! Non orographic Gravity Waves: AF24: removed
1054	!---------------------------------------------
1055
1056	!-----------------------------------------------
1057	! V. Carbon dioxide condensation-sublimation :
1058	!-----------------------------------------------
1059
1060	if (n2cond) then
1061
1062	if (.not.tracer) then
1063	print*,'We need a N2 ice tracer to condense N2'
1064	call abort
1065	endif
1066	call condense_n2(ngrid,nlayer,nq,ptimestep, &
1067	capcal,pplay,pplev,tsurf,pt, &
1068	pdt,zdtsurf,pq,pdq, &
1069	qsurf,zdqsurfc,albedo,emis, &
1070	albedo_bareground,albedo_n2_ice_SPECTV, &
1071	zdtc,zdtsurfc,pdpsrf,zdqc)
1072
1073	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1074	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1075
1076	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1077	dqsurf(1:ngrid,igcm_n2_ice) = dqsurf(1:ngrid,igcm_n2_ice) + zdqsurfc(1:ngrid)
1078
1079	!! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
1080	!! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
1081
1082	! test energy conservation
1083	if(enertest)then
1084	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1085	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
1086	if (is_master) then
1087	print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2'
1088	print*,'In n2cloud surface energy change =',dEtots,' W m-2'
1089	endif
1090	endif
1091
1092	endif ! end of 'n2cond'
1093
1094
1095	!---------------------------------------------
1096	! VI. Specific parameterizations for tracers
1097	!---------------------------------------------
1098
1099	if (tracer) then
1100
1101	! ---------------------
1102	! VI.1. Water and ice !AF24: removed
1103	! ---------------------
1104	! -------------------------
1105	! VI.2. Photochemistry !AF24: removed
1106	! -------------------------
1107	! -------------------------
1108	! VI.3. Aerosol particles
1109	! -------------------------
1110
1111	!Generic Condensation
1112	if (generic_condensation) then
1113	call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, &
1114	pt,pq,pdt,pdq,dt_generic_condensation, &
1115	dqvaplscale_generic,dqcldlscale_generic,rneb_generic)
1116	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dt_generic_condensation(1:ngrid,1:nlayer)
1117	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq)
1118	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq)
1119
1120	if(enertest)then
1121	do ig=1,ngrid
1122	genericconddE(ig) = cppSUM(mass(:,:)dt_generic_condensation(:,:))
1123	enddo
1124
1125	call planetwide_sumval(cppmassarea(:,:)dt_generic_condensation(:,:)/totarea_planet,dEtot)
1126
1127	if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2'
1128	end if
1129
1130	! if (.not. water) then
1131	! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
1132	! Water is the priority
1133	! If you have set water and generic_condensation, then cloudfrac will be water cloudfrac
1134	!
1135	! If you have set generic_condensation (and not water) and you have set several GCS
1136	! then cloudfrac will be only the cloudfrac of the last generic tracer
1137	! (Because it is rewritten every tracer in the loop)
1138	!
1139	! Maybe one should create a cloudfrac_generic(ngrid,nlayer,nq) with 3 dimensions, the last one for tracers
1140
1141	! Let's loop on tracers
1142	cloudfrac(:,:)=0.0
1143	do iq=1,nq
1144	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
1145	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
1146	do l = 1, nlayer
1147	do ig=1,ngrid
1148	cloudfrac(ig,l)=rneb_generic(ig,l,iq)
1149	enddo
1150	enddo
1151	endif
1152	end do ! do iq=1,nq loop on tracers
1153	! endif ! .not. water
1154
1155	endif !generic_condensation
1156
1157	!Generic Rain !AF24: removed
1158
1159	! Sedimentation.
1160	if (sedimentation) then
1161
1162	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1163	zdqssed(1:ngrid,1:nq) = 0.0
1164
1165	! if(watertest)then !AF24: removed
1166
1167	call callsedim(ngrid,nlayer,ptimestep, &
1168	pplev,zzlev,pt,pdt,pq,pdq, &
1169	zdqsed,zdqssed,nq)
1170
1171	! if(watertest)then !AF24: removed
1172
1173	! Whether it falls as rain or snow depends only on the surface temperature
1174	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1175	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1176
1177	!! call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
1178
1179	! ! Test water conservation !AF24: removed
1180
1181	end if ! end of 'sedimentation'
1182
1183	!! call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
1184	!! call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
1185
1186	! ---------------
1187	! VI.4. Updates
1188	! ---------------
1189
1190	! Updating Atmospheric Mass and Tracers budgets.
1191	if(mass_redistrib) then
1192
1193	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
1194	! ( zdqevap(1:ngrid,1:nlayer) &
1195	! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) &
1196	! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) &
1197	! + dqvaplscale(1:ngrid,1:nlayer) )
1198
1199	do ig = 1, ngrid
1200	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1201	enddo
1202
1203	! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1204	! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1205	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1206
1207	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1208	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1209	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1210	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1211
1212	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1213	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1214	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1215	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1216	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1217	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1218	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1219
1220	endif
1221
1222	! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
1223
1224	!! call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
1225	!! call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
1226
1227
1228	! ------------------
1229	! VI.5. Slab Ocean !AF24: removed
1230	! ------------------
1231
1232	! -----------------------------
1233	! VI.6. Surface Tracer Update
1234	! -----------------------------
1235
1236	! AF24: deleted hydrology
1237
1238	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1239
1240	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
1241	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
1242	qsurf_hist(:,:) = qsurf(:,:)
1243
1244	! if(ice_update)then
1245	! ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
1246	! endif
1247
1248	endif! end of if 'tracer'
1249
1250
1251	!------------------------------------------------
1252	! VII. Surface and sub-surface soil temperature
1253	!------------------------------------------------
1254
1255
1256	! ! Increment surface temperature
1257	! if(ok_slab_ocean)then !AF24: removed
1258
1259	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1260
1261	! Compute soil temperatures and subsurface heat flux.
1262	if (callsoil) then
1263	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
1264	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1265	endif
1266
1267
1268	! if (ok_slab_ocean) then !AF24: removed
1269
1270	! Test energy conservation
1271	if(enertest)then
1272	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1273	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1274	endif
1275
1276
1277	!---------------------------------------------------
1278	! VIII. Perform diagnostics and write output files
1279	!---------------------------------------------------
1280
1281	! Note : For output only: the actual model integration is performed in the dynamics.
1282
1283
1284	! Temperature, zonal and meridional winds.
1285	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1286	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1287	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1288
1289	!! Recast thermal plume vertical velocity array for outputs
1290	!! AF24: removed
1291
1292	! Diagnostic.
1293	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1294	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1295
1296	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1297	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1298
1299	if(firstcall)then
1300	zdtdyn(1:ngrid,1:nlayer)=0.0
1301	zdudyn(1:ngrid,1:nlayer)=0.0
1302	endif
1303
1304	! Dynamical heating diagnostic.
1305	do ig=1,ngrid
1306	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1307	enddo
1308
1309	! Tracers.
1310	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1311
1312	! Surface pressure.
1313	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1314
1315
1316
1317	! Surface and soil temperature information
1318	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1319	call planetwide_minval(tsurf(:),Ts2)
1320	call planetwide_maxval(tsurf(:),Ts3)
1321	if(callsoil)then
1322	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1323	if (is_master) then
1324	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1325	print*,Ts1,Ts2,Ts3,TsS
1326	end if
1327	else
1328	if (is_master) then
1329	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1330	print*,Ts1,Ts2,Ts3
1331	endif
1332	end if
1333
1334
1335	! Check the energy balance of the simulation during the run
1336	if(corrk)then
1337
1338	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1339	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1340	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1341	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1342	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1343	do ig=1,ngrid
1344	if(fluxtop_dn(ig).lt.0.0)then
1345	print*,'fluxtop_dn has gone crazy'
1346	print*,'fluxtop_dn=',fluxtop_dn(ig)
1347	print*,'tau_col=',tau_col(ig)
1348	print*,'aerosol=',aerosol(ig,:,:)
1349	print*,'temp= ',pt(ig,:)
1350	print*,'pplay= ',pplay(ig,:)
1351	call abort
1352	endif
1353	end do
1354
1355	if(ngrid.eq.1)then
1356	DYN=0.0
1357	endif
1358
1359	if (is_master) then
1360	print*,' ISR ASR OLR GND DYN [W m^-2]'
1361	print*, ISR,ASR,OLR,GND,DYN
1362	endif
1363
1364	if(enertest .and. is_master)then
1365	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1366	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1367	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1368	endif
1369
1370	if(meanOLR .and. is_master)then
1371	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1372	! to record global radiative balance
1373	open(92,file="rad_bal.out",form='formatted',position='append')
1374	write(92,*) zday,ISR,ASR,OLR
1375	close(92)
1376	open(93,file="tem_bal.out",form='formatted',position='append')
1377	if(callsoil)then
1378	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1379	else
1380	write(93,*) zday,Ts1,Ts2,Ts3
1381	endif
1382	close(93)
1383	endif
1384	endif
1385
1386	endif ! end of 'corrk'
1387
1388
1389	! Diagnostic to test radiative-convective timescales in code.
1390	if(testradtimes)then
1391	open(38,file="tau_phys.out",form='formatted',position='append')
1392	ig=1
1393	do l=1,nlayer
1394	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1395	enddo
1396	close(38)
1397	print*,'As testradtimes enabled,'
1398	print*,'exiting physics on first call'
1399	call abort
1400	endif
1401
1402
1403	! Compute column amounts (kg m-2) if tracers are enabled.
1404	if(tracer)then
1405	qcol(1:ngrid,1:nq)=0.0
1406	do iq=1,nq
1407	do ig=1,ngrid
1408	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1409	enddo
1410	enddo
1411
1412	! Generalised for arbitrary aerosols now. By LK
1413	reffcol(1:ngrid,1:naerkind)=0.0
1414	if(n2cond.and.(iaero_n2.ne.0))then
1415	call n2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_n2))
1416	do ig=1,ngrid
1417	reffcol(ig,iaero_n2) = SUM(zq(ig,1:nlayer,igcm_n2_ice)reffrad(ig,1:nlayer,iaero_n2)mass(ig,1:nlayer))
1418	enddo
1419	endif
1420	! if(water.and.(iaero_h2o.ne.0))then !AF24: removed
1421
1422	endif ! end of 'tracer'
1423
1424
1425	! ! Test for water conservation. !AF24: removed
1426
1427	! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
1428	if(generic_condensation)then
1429	RH_generic(:,:,:)=0.0
1430	do iq=1,nq
1431
1432	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
1433
1434	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
1435
1436	do l = 1, nlayer
1437	do ig=1,ngrid
1438	call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
1439	RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_vap) / qsat_generic(ig,l,iq)
1440	enddo
1441	enddo
1442
1443	end if
1444
1445	end do ! iq=1,nq
1446
1447	endif ! end of 'generic_condensation'
1448
1449
1450	if (is_master) print,'--> Ls =',zls180./pi
1451
1452
1453	!----------------------------------------------------------------------
1454	! Writing NetCDF file "RESTARTFI" at the end of the run
1455	!----------------------------------------------------------------------
1456
1457	! Note: 'restartfi' is stored just before dynamics are stored
1458	! in 'restart'. Between now and the writting of 'restart',
1459	! there will have been the itau=itau+1 instruction and
1460	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1461	! thus we store for time=time+dtvr
1462
1463
1464
1465	if(lastcall) then
1466	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1467
1468	!! Update surface ice distribution to iterate to steady state if requested
1469	!! AF24: removed
1470
1471	! endif
1472	#ifndef MESOSCALE
1473
1474	if (ngrid.ne.1) then
1475	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1476
1477	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
1478	ptimestep,ztime_fin, &
1479	tsurf,tsoil,emis,q2,qsurf_hist)
1480	endif
1481	#endif
1482	! if(ok_slab_ocean) then
1483	! call ocean_slab_final!(tslab, seaice)
1484	! end if
1485
1486	endif ! end of 'lastcall'
1487
1488
1489	! -----------------------------------------------------------------
1490	! WSTATS: Saving statistics
1491	! -----------------------------------------------------------------
1492	! ("stats" stores and accumulates key variables in file "stats.nc"
1493	! which can later be used to make the statistic files of the run:
1494	! if flag "callstats" from callphys.def is .true.)
1495
1496
1497	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1498	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1499	call wstats(ngrid,"fluxsurf_lw", &
1500	"Thermal IR radiative flux to surface","W.m-2",2, &
1501	fluxsurf_lw)
1502	call wstats(ngrid,"fluxtop_lw", &
1503	"Thermal IR radiative flux to space","W.m-2",2, &
1504	fluxtop_lw)
1505
1506	! call wstats(ngrid,"fluxsurf_sw", &
1507	! "Solar radiative flux to surface","W.m-2",2, &
1508	! fluxsurf_sw_tot)
1509	! call wstats(ngrid,"fluxtop_sw", &
1510	! "Solar radiative flux to space","W.m-2",2, &
1511	! fluxtop_sw_tot)
1512
1513
1514	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1515	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1516	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1517	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
1518	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
1519	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
1520	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1521	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1522	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1523	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1524	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1525
1526	if (tracer) then
1527	do iq=1,nq
1528	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1529	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1530	'kg m^-2',2,qsurf(1,iq) )
1531	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
1532	'kg m^-2',2,qcol(1,iq) )
1533
1534	! call wstats(ngrid,trim(noms(iq))//'_reff', &
1535	! trim(noms(iq))//'_reff', &
1536	! 'm',3,reffrad(1,1,iq))
1537
1538	end do
1539
1540	endif ! end of 'tracer'
1541
1542	!AF24: deleted slab ocean and water
1543
1544	if(lastcall.and.callstats) then
1545	write (,) "Writing stats..."
1546	call mkstats(ierr)
1547	endif
1548
1549
1550	#ifndef MESOSCALE
1551
1552	!-----------------------------------------------------------------------------------------------------
1553	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
1554	!
1555	! Note 1 : output with period "ecritphy", set in "run.def"
1556	!
1557	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
1558	! but its preferable to keep all the calls in one place ...
1559	!-----------------------------------------------------------------------------------------------------
1560
1561	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
1562	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
1563	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
1564	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
1565	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1566	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1567	call writediagfi(ngrid,"temp","temperature","K",3,zt)
1568	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
1569	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1570	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1571	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1572	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
1573
1574	! Subsurface temperatures
1575	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1576	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
1577
1578	! ! Oceanic layers !AF24: removed
1579
1580	! ! Thermal plume model !AF24: removed
1581
1582	! GW non-oro outputs !AF24: removed
1583
1584	! Total energy balance diagnostics
1585	if(callrad)then
1586
1587	!call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
1588	!call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
1589	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1590	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1591	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1592	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
1593
1594	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
1595	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
1596	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
1597	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
1598	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
1599	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
1600
1601	! if(ok_slab_ocean) then
1602	! call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
1603	! else
1604	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
1605	! endif
1606
1607	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1608
1609	endif ! end of 'callrad'
1610
1611	if(enertest) then
1612
1613	if (calldifv) then
1614
1615	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
1616	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
1617
1618	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
1619	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
1620	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
1621
1622	endif
1623
1624	if (corrk) then
1625	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
1626	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
1627	endif
1628
1629	! if(watercond) then !AF24: removed
1630
1631	if (generic_condensation) then
1632
1633	call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
1634	call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
1635
1636	endif
1637
1638	endif ! end of 'enertest'
1639
1640	! Diagnostics of optical thickness
1641	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
1642	if (diagdtau) then
1643	do nw=1,L_NSPECTV
1644	write(str2,'(i2.2)') nw
1645	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
1646	enddo
1647	do nw=1,L_NSPECTI
1648	write(str2,'(i2.2)') nw
1649	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
1650	enddo
1651	endif
1652
1653
1654	! Temporary inclusions for heating diagnostics.
1655	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1656	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1657	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
1658	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1659
1660	! For Debugging.
1661	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
1662	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
1663
1664
1665	! Output aerosols.
1666	if (igcm_n2_ice.ne.0.and.iaero_n2.ne.0) &
1667	call writediagfi(ngrid,'N2ice_reff','N2ice_reff','m',3,reffrad(1,1,iaero_n2))
1668	if (igcm_n2_ice.ne.0.and.iaero_n2.ne.0) &
1669	call writediagfi(ngrid,'N2ice_reffcol','N2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_n2))
1670	! if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & !AF24: removed
1671
1672	! Output tracers.
1673	if (tracer) then
1674
1675	do iq=1,nq
1676	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1677	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1678	'kg m^-2',2,qsurf_hist(1,iq) )
1679	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
1680	'kg m^-2',2,qcol(1,iq) )
1681	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1682	! 'kg m^-2',2,qsurf(1,iq) )
1683
1684	! if(watercond.or.CLFvarying)then !AF24: removed
1685
1686	if(generic_condensation)then
1687	call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
1688	call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
1689	call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
1690	endif
1691
1692	! if(generic_rain)then !AF24: removed
1693	! if((hydrology).and.(.not.ok_slab_ocean))then !AF24: removed
1694
1695	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
1696
1697	enddo ! end of 'nq' loop
1698
1699	endif ! end of 'tracer'
1700
1701
1702	! Output spectrum.
1703	if(specOLR.and.corrk)then
1704	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
1705	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
1706	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
1707	endif
1708
1709	#else
1710	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
1711	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
1712	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
1713	if (.not.calldifv) comm_LATENT_HF(:)=0.0
1714	! if ((tracer).and.(water)) then !AF24: removed
1715
1716	if ((tracer).and.(generic_condensation)) then
1717	! .and.(.not. water)
1718
1719	! If you have set generic_condensation (and not water) and you have set several GCS
1720	! then the outputs given to WRF will be only the ones for the last generic tracer
1721	! (Because it is rewritten every tracer in the loop)
1722	! WRF can take only one moist tracer
1723
1724	do iq=1,nq
1725	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
1726
1727	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
1728
1729	reffrad_generic_zeros_for_wrf(:,:) = 1.
1730
1731	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
1732	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
1733	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
1734	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_vap)
1735	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
1736	! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
1737	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
1738	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
1739	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
1740	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
1741	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
1742
1743	endif
1744	end do ! do iq=1,nq loop on tracers
1745
1746	else
1747	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
1748	comm_TOTCLOUDFRAC(1:ngrid)=0.
1749	comm_SURFRAIN(1:ngrid)=0.
1750	comm_DQVAP(1:ngrid,1:nlayer)=0.
1751	comm_DQICE(1:ngrid,1:nlayer)=0.
1752	! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
1753	comm_REEVAP(1:ngrid)=0.
1754	comm_DTRAIN(1:ngrid,1:nlayer)=0.
1755	comm_DTLSC(1:ngrid,1:nlayer)=0.
1756	comm_RH(1:ngrid,1:nlayer)=0.
1757
1758	endif ! if water, if generic_condensation, else
1759
1760	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
1761	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
1762	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
1763	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
1764	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
1765	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
1766	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
1767	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
1768
1769	! if (turb_resolved) then
1770	! open(17,file='lsf.txt',form='formatted',status='old')
1771	! rewind(17)
1772	! DO l=1,nlayer
1773	! read(17,*) lsf_dt(l),lsf_dq(l)
1774	! ENDDO
1775	! close(17)
1776	! do ig=1,ngrid
1777	! if ((tracer).and.(water)) then
1778	! pdq(ig,:,igcm_h2o_vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:)
1779	! endif
1780	! pdt(ig,:) = pdt(ig,:) + lsf_dt(:)
1781	! comm_HR_DYN(ig,:) = lsf_dt(:)
1782	! enddo
1783	! endif
1784	#endif
1785
1786	! XIOS outputs
1787	#ifdef CPP_XIOS
1788	! Send fields to XIOS: (NB these fields must also be defined as
1789	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1790	CALL send_xios_field("ls",zls)
1791
1792	CALL send_xios_field("ps",ps)
1793	CALL send_xios_field("area",cell_area)
1794	CALL send_xios_field("p",pplay)
1795	CALL send_xios_field("temperature",zt)
1796	CALL send_xios_field("u",zu)
1797	CALL send_xios_field("v",zv)
1798	CALL send_xios_field("omega",omega)
1799
1800	! IF (calltherm) THEN !AF24: removed
1801	! IF (water) THEN !AF24: removed
1802
1803	CALL send_xios_field("ISR",fluxtop_dn)
1804	CALL send_xios_field("OLR",fluxtop_lw)
1805	CALL send_xios_field("ASR",fluxabs_sw)
1806
1807	if (specOLR .and. corrk) then
1808	call send_xios_field("OLR3D",OLR_nu)
1809	call send_xios_field("IR_Bandwidth",DWNI)
1810	call send_xios_field("VI_Bandwidth",DWNV)
1811	call send_xios_field("OSR3D",OSR_nu)
1812	call send_xios_field("GSR3D",GSR_nu)
1813	endif
1814
1815	if (lastcall.and.is_omp_master) then
1816	write(,) "physiq: call xios_context_finalize"
1817	call xios_context_finalize
1818	endif
1819	#endif
1820
1821	if (check_physics_outputs) then
1822	! Check the validity of updated fields at the end of the physics step
1823	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
1824	endif
1825
1826	icount=icount+1
1827
1828	end subroutine physiq
1829
1830	end module physiq_mod

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

Download in other formats:

Original Format