Context Navigation

physiq.F90 @ 1477

Last change on this file since 1477 was 1477, checked in by mturbet, 9 years ago
Cleaning of physiq.F90. Condense_cloud becomes Condense_co2
Property svn:executable set to ``*
File size: 86.6 KB

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

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

Download in other formats:

Original Format