Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90 @ 227

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

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

Download in other formats: