Context Navigation

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

Last change on this file since 232 was 232, checked in by ymipsl, 10 years ago

Work around for init etat0 state => missing field in the startfi are initialize to constant value.

YM

Property svn:executable set to *

File size: 94.2 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	!ym call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
562	!ym day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
563	!ym cloudfrac,totcloudfrac,hice, &
564	!ym rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
565
566	call phyetat0_academic(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
567	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
568	cloudfrac,totcloudfrac,hice, &
569	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
570
571	if (pday.ne.day_ini) then
572	write(,) "ERROR in physiq.F90:"
573	write(,) "bad synchronization between physics and dynamics"
574	write(,) "dynamics day: ",pday
575	write(,) "physics day: ",day_ini
576	write(,) "taking dynamics day for physics: ",day_ini
577	day_ini=pday
578	!ym stop
579	endif
580
581	write (,) 'In physiq day_ini =', day_ini
582
583	! Initialize albedo and orbital calculation
584	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
585	call surfini(ngrid,nq,qsurf,albedo0)
586	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
587
588	albedo(:)=albedo0(:)
589
590	if(tlocked)then
591	print*,'Planet is tidally locked at resonance n=',nres
592	print*,'Make sure you have the right rotation rate!!!'
593	endif
594
595	! Initialize oceanic variables
596	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
597
598	if (ok_slab_ocean)then
599
600	call ocean_slab_init(ngrid,ptimestep, tslab, &
601	sea_ice, pctsrf_sic)
602
603	knindex(:) = 0
604
605	do ig=1,ngrid
606	zmasq(ig)=1
607	knindex(ig) = ig
608	if (nint(rnat(ig)).eq.0) then
609	zmasq(ig)=0
610	endif
611	enddo
612
613
614	CALL init_masquv(ngrid,zmasq)
615
616
617	endif
618
619
620	! initialize soil
621	! ~~~~~~~~~~~~~~~
622	if (callsoil) then
623	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
624	ptimestep,tsurf,tsoil,capcal,fluxgrd)
625
626	if (ok_slab_ocean) then
627	do ig=1,ngrid
628	if (nint(rnat(ig)).eq.2) then
629	capcal(ig)=capcalocean
630	if (pctsrf_sic(ig).gt.0.5) then
631	capcal(ig)=capcalseaice
632	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
633	capcal(ig)=capcalsno
634	endif
635	endif
636	endif
637	enddo
638	endif
639
640
641
642
643	else
644	print*,'WARNING! Thermal conduction in the soil turned off'
645	capcal(:)=1.e6
646	fluxgrd(:)=intheat
647	print*,'Flux from ground = ',intheat,' W m^-2'
648	endif
649	icount=1
650
651	! decide whether to update ice at end of run
652	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
653	ice_update=.false.
654	if(sourceevol)then
655	!$OMP MASTER
656	open(128,file='num_run',form='formatted', &
657	status="old",iostat=ierr)
658	if (ierr.ne.0) then
659	write(,) "physiq: Error! No num_run file!"
660	write(,) " (which is needed for sourceevol option)"
661	stop
662	endif
663	read(128,*) num_run
664	close(128)
665	!$OMP END MASTER
666	!$OMP BARRIER
667
668	if(num_run.ne.0.and.mod(num_run,2).eq.0)then
669	!if(num_run.ne.0.and.mod(num_run,3).eq.0)then
670	print*,'Updating ice at end of this year!'
671	ice_update=.true.
672	ice_min(:)=1.e4
673	endif
674	endif
675
676	! In newstart now, will have to be remove (BC 2014)
677	! define surface as continent or ocean
678	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
679	if (.not.ok_slab_ocean) then
680	rnat(:)=1.
681	do ig=1,ngrid
682	! if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
683	if(inertiedat(ig,1).gt.1.E4)then
684	rnat(ig)=0
685	endif
686	enddo
687
688	print*,'WARNING! Surface type currently decided by surface inertia'
689	print*,'This should be improved e.g. in newstart.F'
690	endif!(.not.ok_slab_ocean)
691
692	! initialise surface history variable
693	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
694	qsurf_hist(:,:)=qsurf(:,:)
695
696	! initialise variable for dynamical heating diagnostic
697	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
698	ztprevious(:,:)=pt(:,:)
699
700	! Set temperature just above condensation temperature (for Early Mars)
701	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
702	if(nearco2cond) then
703	write(,)' WARNING! Starting at Tcond+1K'
704	do l=1, nlayer
705	do ig=1,ngrid
706	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
707	-pt(ig,l)) / ptimestep
708	enddo
709	enddo
710	endif
711
712	if(meanOLR)then
713	! to record global radiative balance
714	call system('rm -f rad_bal.out')
715	! to record global mean/max/min temperatures
716	call system('rm -f tem_bal.out')
717	! to record global hydrological balance
718	call system('rm -f h2o_bal.out')
719	endif
720
721	watertest=.false.
722	if(water)then
723	! initialise variables for water cycle
724
725	if(enertest)then
726	watertest = .true.
727	endif
728
729	if(ice_update)then
730	ice_initial(:)=qsurf(:,igcm_h2o_ice)
731	endif
732
733	endif
734	call su_watercycle ! even if we don't have a water cycle, we might
735	! need epsi for the wvp definitions in callcorrk.F
736
737	if (ngrid.ne.1) then ! no need to create a restart file in 1d
738	call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
739	ptimestep,pday+nday,time_phys,area, &
740	albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
741	endif
742
743	endif ! (end of "if firstcall")
744
745	! ---------------------------------------------------
746	! 1.2 Initializations done at every physical timestep:
747	! ---------------------------------------------------
748	! Initialize various variables
749	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
750
751	pdu(1:ngrid,1:nlayer) = 0.0
752	pdv(1:ngrid,1:nlayer) = 0.0
753	if ( .not.nearco2cond ) then
754	pdt(1:ngrid,1:nlayer) = 0.0
755	endif
756	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
757	pdpsrf(1:ngrid) = 0.0
758	zflubid(1:ngrid) = 0.0
759	zdtsurf(1:ngrid) = 0.0
760	dqsurf(1:ngrid,1:nq)= 0.0
761	flux_sens_lat(1:ngrid) = 0.0
762	taux(1:ngrid) = 0.0
763	tauy(1:ngrid) = 0.0
764
765
766
767
768	zday=pday+ptime ! compute time, in sols (and fraction thereof)
769
770	! Compute Stellar Longitude (Ls)
771	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
772	if (season) then
773	call stellarlong(zday,zls)
774	else
775	call stellarlong(float(day_ini),zls)
776	end if
777
778
779
780	! Compute variations of g with latitude (oblate case)
781	!
782	if (oblate .eqv. .false.) then
783	glat(:) = g
784	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
785	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
786
787	call abort
788	else
789
790	gmplanet = MassPlanetgrav1e24
791	do ig=1,ngrid
792	glat(ig)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - lati(ig))))
793	end do
794	endif
795
796	!! write(,) 'lati, glat =',lati(1)180./pi,glat(1), lati(ngrid/2)180./pi, glat(ngrid/2)
797
798
799
800	! Compute geopotential between layers
801	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
802
803	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
804	do l=1,nlayer
805	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
806	enddo
807
808	zzlev(1:ngrid,1)=0.
809	zzlev(1:ngrid,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
810
811	do l=2,nlayer
812	do ig=1,ngrid
813	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
814	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
815	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
816	enddo
817	enddo
818	! Potential temperature calculation may not be the same in physiq and dynamic...
819
820	! Compute potential temperature
821	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
822
823	do l=1,nlayer
824	do ig=1,ngrid
825	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
826	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
827	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
828	massarea(ig,l)=mass(ig,l)*area(ig)
829	enddo
830	enddo
831
832	! Compute vertical velocity (m/s) from vertical mass flux
833	! w = F / (rhoarea) and rho = rT/P
834	! but first linearly interpolate mass flux to mid-layers
835	do l=1,nlayer-1
836	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
837	enddo
838	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
839	do l=1,nlayer
840	pw(1:ngrid,l)=(pw(1:ngrid,l)*pplay(1:ngrid,l)) / &
841	(rpt(1:ngrid,l)area(1:ngrid))
842	enddo
843
844	!-----------------------------------------------------------------------
845	! 2. Compute radiative tendencies
846	!-----------------------------------------------------------------------
847
848	if (callrad) then
849	if( mod(icount-1,iradia).eq.0.or.lastcall) then
850
851	! Compute local stellar zenith angles
852	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
853	call orbite(zls,dist_star,declin)
854
855	if (tlocked) then
856	ztim1=SIN(declin)
857	! ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
858	! ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
859	! JL12 corrects tidally resonant cases. nres=omega_rot/omega_orb
860	ztim2=COS(declin)COS(2.pi(zday/year_day)nres - zls)
861	ztim3=-COS(declin)SIN(2.pi(zday/year_day)nres - zls)
862
863	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
864	ztim1,ztim2,ztim3,mu0,fract, flatten)
865
866	elseif (diurnal) then
867	ztim1=SIN(declin)
868	ztim2=COS(declin)COS(2.pi*(zday-.5))
869	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
870
871	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
872	ztim1,ztim2,ztim3,mu0,fract, flatten)
873	else if(diurnal .eqv. .false.) then
874
875	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad,flatten)
876	! WARNING: this function appears not to work in 1D
877
878	endif
879
880	!! Eclipse incoming sunlight (e.g. Saturn ring shadowing)
881
882	if(rings_shadow) then
883	write(,) 'Rings shadow activated'
884	if(diurnal .eqv. .false.) then ! we need to compute the daily average insolation
885	pas = 1./nb_hours
886	ptime_day = 0.
887	fract(:) = 0.
888	ALLOCATE(tmp_fract(ngrid))
889	ALLOCATE(tmp_mu0(ngrid))
890	tmp_fract(:) = 0.
891	eclipse(:) = 0.
892	tmp_mu0(:) = 0.
893
894	do m=1, nb_hours
895	ptime_day = m*pas
896	call stellarlong(pday+ptime_day,tmp_zls)
897	call orbite(tmp_zls,dist_star,declin)
898	ztim1=SIN(declin)
899	ztim2=COS(declin)COS(2.pi*(pday+ptime_day-.5))
900	ztim3=-COS(declin)SIN(2.pi*(pday+ptime_day-.5))
901
902	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
903	ztim1,ztim2,ztim3,tmp_mu0,tmp_fract, flatten)
904
905	call rings(ngrid, declin, ptime_day, rad, flatten, eclipse)
906
907	fract(:) = fract(:) + (1.-eclipse(:))*tmp_fract(:) !! fract prend en compte l'ombre des anneaux et l'alternance jour/nuit
908	enddo
909
910	DEALLOCATE(tmp_fract)
911	DEALLOCATE(tmp_mu0)
912
913	fract(:) = fract(:)/nb_hours
914
915	else
916	call rings(ngrid, declin, ptime, rad, 0., eclipse)
917	fract(:) = fract(:) * (1.-eclipse)
918	endif
919	else
920	eclipse(:) = 0.0
921	endif
922
923	if (corrk) then
924
925	! a) Call correlated-k radiative transfer scheme
926	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
927
928	if(kastprof)then
929	print*,'kastprof should not = true here'
930	call abort
931	endif
932	if(water) then
933	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
934	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
935	! take into account water effect on mean molecular weight
936	else
937	muvar(1:ngrid,1:nlayer+1)=mugaz
938	endif
939
940	! standard callcorrk
941
942
943
944
945	if(ok_slab_ocean) then
946	tsurf2(:)=tsurf(:)
947	do ig=1,ngrid
948	if (nint(rnat(ig))==0) then
949	tsurf(ig)=((1.-pctsrf_sic(ig))tslab(ig,1)4+pctsrf_sic(ig)tsea_ice(ig)4)0.25
950	endif
951	enddo
952	endif!(ok_slab_ocean)
953
954
955	clearsky=.false.
956	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
957	albedo,emis,mu0,pplev,pplay,pt, &
958	tsurf,fract,dist_star,aerosol,muvar, &
959	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
960	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
961	tau_col,cloudfrac,totcloudfrac, &
962	clearsky,firstcall,lastcall)
963
964	! Option to call scheme once more for clear regions
965	if(CLFvarying)then
966
967	! ---> PROBLEMS WITH ALLOCATED ARRAYS
968	! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
969	clearsky=.true.
970	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
971	albedo,emis,mu0,pplev,pplay,pt, &
972	tsurf,fract,dist_star,aerosol,muvar, &
973	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
974	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, &
975	tau_col1,cloudfrac,totcloudfrac, &
976	clearsky,firstcall,lastcall)
977	clearsky = .false. ! just in case.
978
979
980	! Sum the fluxes and heating rates from cloudy/clear cases
981	do ig=1,ngrid
982	tf=totcloudfrac(ig)
983	ntf=1.-tf
984
985	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
986	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
987	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
988	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
989	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
990
991	zdtlw(ig,1:nlayer) = ntfzdtlw1(ig,1:nlayer) + tfzdtlw(ig,1:nlayer)
992	zdtsw(ig,1:nlayer) = ntfzdtsw1(ig,1:nlayer) + tfzdtsw(ig,1:nlayer)
993
994	OSR_nu(ig,1:L_NSPECTV) = ntfOSR_nu1(ig,1:L_NSPECTV) + tfOSR_nu(ig,1:L_NSPECTV)
995	OLR_nu(ig,1:L_NSPECTI) = ntfOLR_nu1(ig,1:L_NSPECTI) + tfOLR_nu(ig,1:L_NSPECTI)
996
997	enddo
998
999	endif !CLFvarying
1000
1001	if(ok_slab_ocean) then
1002	tsurf(:)=tsurf2(:)
1003	endif!(ok_slab_ocean)
1004
1005
1006	! Radiative flux from the sky absorbed by the surface (W.m-2)
1007	GSR=0.0
1008	fluxrad_sky(1:ngrid)=emis(1:ngrid)fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)(1.-albedo(1:ngrid))
1009
1010	!if(noradsurf)then ! no lower surface; SW flux just disappears
1011	! GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
1012	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1013	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1014	!endif
1015
1016	! Net atmospheric radiative heating rate (K.s-1)
1017	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1018
1019	elseif(newtonian)then
1020
1021	! b) Call Newtonian cooling scheme
1022	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1023	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
1024
1025	zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
1026	! e.g. surface becomes proxy for 1st atmospheric layer ?
1027
1028	else
1029
1030	! c) Atmosphere has no radiative effect
1031	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1032	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
1033	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1034	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
1035	endif
1036	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1037	fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
1038	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
1039	! radiation skips the atmosphere entirely
1040
1041
1042	dtrad(1:ngrid,1:nlayer)=0.0
1043	! hence no atmospheric radiative heating
1044
1045	endif ! if corrk
1046
1047	endif ! of if(mod(icount-1,iradia).eq.0)
1048
1049
1050	! Transformation of the radiative tendencies
1051	! ------------------------------------------
1052
1053	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1054	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
1055	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1056	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1057
1058	!-------------------------
1059	! test energy conservation
1060	if(enertest)then
1061	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
1062	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
1063	call planetwide_sumval(fluxsurf_sw(:)(1.-albedo(:))area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1064	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))area(:)/totarea_planet,dEtotsLW)
1065	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
1066	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
1067	if (is_master) then
1068	print*,'---------------------------------------------------------------'
1069	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
1070	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
1071	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
1072	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
1073	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
1074	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
1075	endif
1076	endif
1077	!-------------------------
1078
1079	endif ! of if (callrad)
1080
1081	!-----------------------------------------------------------------------
1082	! 4. Vertical diffusion (turbulent mixing):
1083	! -----------------------------------------
1084
1085	if (calldifv) then
1086
1087	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1088
1089	zdum1(1:ngrid,1:nlayer)=0.0
1090	zdum2(1:ngrid,1:nlayer)=0.0
1091
1092
1093	!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff
1094	if (UseTurbDiff) then
1095
1096	call turbdiff(ngrid,nlayer,nq,rnat, &
1097	ptimestep,capcal,lwrite, &
1098	pplay,pplev,zzlay,zzlev,z0, &
1099	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
1100	zdum1,zdum2,pdt,pdq,zflubid, &
1101	zdudif,zdvdif,zdtdif,zdtsdif, &
1102	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
1103	taux,tauy,lastcall)
1104
1105	else
1106
1107	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1108
1109	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
1110	ptimestep,capcal,lwrite, &
1111	pplay,pplev,zzlay,zzlev,z0, &
1112	pu,pv,zh,pq,tsurf,emis,qsurf, &
1113	zdum1,zdum2,zdh,pdq,zflubid, &
1114	zdudif,zdvdif,zdhdif,zdtsdif, &
1115	sensibFlux,q2,zdqdif,zdqsdif, &
1116	taux,tauy,lastcall)
1117
1118	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1119	zdqevap(1:ngrid,1:nlayer)=0.
1120
1121	end if !turbdiff
1122
1123	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1124	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1125	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1126	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1127
1128	if(ok_slab_ocean)then
1129	flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
1130	endif
1131
1132
1133
1134	if (tracer) then
1135	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1136	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1137	end if ! of if (tracer)
1138
1139	!-------------------------
1140	! test energy conservation
1141	if(enertest)then
1142	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
1143	do ig = 1, ngrid
1144	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1145	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1146	enddo
1147	call planetwide_sumval(dEdiff(:)*area(:)/totarea_planet,dEtot)
1148	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1149	call planetwide_sumval(dEdiffs(:)*area(:)/totarea_planet,dEtots)
1150	call planetwide_sumval(sensibFlux(:)*area(:)/totarea_planet,AtmToSurf_TurbFlux)
1151	if (is_master) then
1152	if (UseTurbDiff) then
1153	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1154	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1155	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1156	else
1157	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1158	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1159	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1160	end if
1161	endif ! of if (is_master)
1162	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
1163	! but not given back elsewhere
1164	endif
1165	!-------------------------
1166
1167	!-------------------------
1168	! test water conservation
1169	if(watertest.and.water)then
1170	call planetwide_sumval(massarea(:,:)zdqdif(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dWtot_tmp)
1171	call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)area(:)ptimestep/totarea_planet,dWtots_tmp)
1172	do ig = 1, ngrid
1173	vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
1174	Enddo
1175	call planetwide_sumval(massarea(:,:)zdqdif(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1176	call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)area(:)ptimestep/totarea_planet,dWtots)
1177	dWtot = dWtot + dWtot_tmp
1178	dWtots = dWtots + dWtots_tmp
1179	do ig = 1, ngrid
1180	vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
1181	Enddo
1182	call planetwide_maxval(vdifcncons(:),nconsMAX)
1183
1184	if (is_master) then
1185	print*,'---------------------------------------------------------------'
1186	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
1187	print*,'In difv surface water change =',dWtots,' kg m-2'
1188	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
1189	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
1190	endif
1191
1192	endif
1193	!-------------------------
1194
1195	else
1196
1197	if(.not.newtonian)then
1198
1199	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1200
1201	endif
1202
1203	endif ! of if (calldifv)
1204
1205
1206	!-----------------------------------------------------------------------
1207	! 5. Dry convective adjustment:
1208	! -----------------------------
1209
1210	if(calladj) then
1211
1212	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1213	zduadj(1:ngrid,1:nlayer)=0.0
1214	zdvadj(1:ngrid,1:nlayer)=0.0
1215	zdhadj(1:ngrid,1:nlayer)=0.0
1216
1217
1218	call convadj(ngrid,nlayer,nq,ptimestep, &
1219	pplay,pplev,zpopsk, &
1220	pu,pv,zh,pq, &
1221	pdu,pdv,zdh,pdq, &
1222	zduadj,zdvadj,zdhadj, &
1223	zdqadj)
1224
1225	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1226	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1227	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1228	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1229
1230	if(tracer) then
1231	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1232	end if
1233
1234	!-------------------------
1235	! test energy conservation
1236	if(enertest)then
1237	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1238	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1239	endif
1240	!-------------------------
1241
1242	!-------------------------
1243	! test water conservation
1244	if(watertest)then
1245	call planetwide_sumval(massarea(:,:)zdqadj(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dWtot_tmp)
1246	do ig = 1, ngrid
1247	cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
1248	Enddo
1249	call planetwide_sumval(massarea(:,:)zdqadj(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1250	dWtot = dWtot + dWtot_tmp
1251	do ig = 1, ngrid
1252	cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
1253	Enddo
1254	call planetwide_maxval(cadjncons(:),nconsMAX)
1255
1256	if (is_master) then
1257	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
1258	print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
1259	endif
1260	endif
1261	!-------------------------
1262
1263	endif ! of if(calladj)
1264
1265	!-----------------------------------------------------------------------
1266	! 6. Carbon dioxide condensation-sublimation:
1267	! -------------------------------------------
1268
1269	if (co2cond) then
1270	if (.not.tracer) then
1271	print*,'We need a CO2 ice tracer to condense CO2'
1272	call abort
1273	endif
1274	call condense_cloud(ngrid,nlayer,nq,ptimestep, &
1275	capcal,pplay,pplev,tsurf,pt, &
1276	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1277	qsurf(1,igcm_co2_ice),albedo,emis, &
1278	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1279	zdqc)
1280
1281	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1282	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
1283	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
1284	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1285
1286	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1287	! Note: we do not add surface co2ice tendency
1288	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
1289
1290	!-------------------------
1291	! test energy conservation
1292	if(enertest)then
1293	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1294	call planetwide_sumval(capcal(:)zdtsurfc(:)area(:)/totarea_planet,dEtots)
1295	if (is_master) then
1296	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
1297	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
1298	endif
1299	endif
1300	!-------------------------
1301
1302	endif ! of if (co2cond)
1303
1304
1305	!-----------------------------------------------------------------------
1306	! 7. Specific parameterizations for tracers
1307	! -----------------------------------------
1308
1309	if (tracer) then
1310
1311	! 7a. Water and ice
1312	! ---------------
1313	if (water) then
1314
1315	! ----------------------------------------
1316	! Water ice condensation in the atmosphere
1317	! ----------------------------------------
1318	if(watercond.and.(RLVTT.gt.1.e-8))then
1319
1320	! ----------------
1321	! Moist convection
1322	! ----------------
1323
1324	dqmoist(1:ngrid,1:nlayer,1:nq)=0.
1325	dtmoist(1:ngrid,1:nlayer)=0.
1326
1327	call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1328
1329	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
1330	+dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
1331	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
1332	+dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
1333	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
1334
1335	!-------------------------
1336	! test energy conservation
1337	if(enertest)then
1338	call planetwide_sumval(cppmassarea(:,:)dtmoist(:,:)/totarea_planet,dEtot)
1339	call planetwide_maxval(dtmoist(:,:),dtmoist_max)
1340	call planetwide_minval(dtmoist(:,:),dtmoist_min)
1341	madjdEz(:,:)=cppmass(:,:)dtmoist(:,:)
1342	do ig=1,ngrid
1343	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
1344	enddo
1345	if (is_master) then
1346	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
1347	print,'In moistadj MAX atmospheric energy change =',dtmoist_maxptimestep,'K/step'
1348	print,'In moistadj MIN atmospheric energy change =',dtmoist_minptimestep,'K/step'
1349	endif
1350
1351	! test energy conservation
1352	call planetwide_sumval(massarea(:,:)dqmoist(:,:,igcm_h2o_vap)ptimestep/totarea_planet+ &
1353	massarea(:,:)dqmoist(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1354	if (is_master) print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
1355	endif
1356	!-------------------------
1357
1358
1359	! --------------------------------
1360	! Large scale condensation/evaporation
1361	! --------------------------------
1362	call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
1363
1364	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer)
1365	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer)
1366	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
1367
1368	!-------------------------
1369	! test energy conservation
1370	if(enertest)then
1371	lscaledEz(:,:) = cppmass(:,:)dtlscale(:,:)
1372	do ig=1,ngrid
1373	lscaledE(ig) = cppSUM(mass(:,:)dtlscale(:,:))
1374	enddo
1375	call planetwide_sumval(cppmassarea(:,:)dtlscale(:,:)/totarea_planet,dEtot)
1376	! if(isnan(dEtot)) then ! NB: isnan() is not a standard function...
1377	! print*,'Nan in largescale, abort'
1378	! STOP
1379	! endif
1380	if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
1381
1382	! test water conservation
1383	call planetwide_sumval(massarea(:,:)dqvaplscale(:,:)ptimestep/totarea_planet+ &
1384	SUM(massarea(:,:)dqcldlscale(:,:))ptimestep/totarea_planet,dWtot)
1385	if (is_master) print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
1386	endif
1387	!-------------------------
1388
1389	! compute cloud fraction
1390	do l = 1, nlayer
1391	do ig=1,ngrid
1392	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1393	enddo
1394	enddo
1395
1396	endif ! of if (watercondense)
1397
1398
1399	! --------------------------------
1400	! Water ice / liquid precipitation
1401	! --------------------------------
1402	if(waterrain)then
1403
1404	zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0
1405	zdqsrain(1:ngrid) = 0.0
1406	zdqssnow(1:ngrid) = 0.0
1407
1408	call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
1409	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
1410
1411	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
1412	+zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
1413	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
1414	+zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
1415	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
1416	dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) ! a bug was here
1417	dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
1418	rainout(1:ngrid) = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic
1419
1420	!-------------------------
1421	! test energy conservation
1422	if(enertest)then
1423	call planetwide_sumval(cppmassarea(:,:)zdtrain(:,:)/totarea_planet,dEtot)
1424	if (is_master) print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
1425	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_ice)/totarea_planetRLVTT/cpp,dItot_tmp)
1426	call planetwide_sumval(area(:)zdqssnow(:)/totarea_planetRLVTT/cpp,dItot)
1427	dItot = dItot + dItot_tmp
1428	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dVtot_tmp)
1429	call planetwide_sumval(area(:)zdqsrain(:)/totarea_planetRLVTT/cpp,dVtot)
1430	dVtot = dVtot + dVtot_tmp
1431	dEtot = dItot + dVtot
1432	if (is_master) then
1433	print*,'In rain dItot =',dItot,' W m-2'
1434	print*,'In rain dVtot =',dVtot,' W m-2'
1435	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
1436	endif
1437
1438	! test water conservation
1439	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_vap)ptimestep/totarea_planet+ &
1440	massarea(:,:)zdqrain(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1441	call planetwide_sumval((zdqsrain(:)+zdqssnow(:))area(:)ptimestep/totarea_planet,dWtots)
1442	if (is_master) then
1443	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
1444	print*,'In rain surface water change =',dWtots,' kg m-2'
1445	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
1446	endif
1447	endif
1448	!-------------------------
1449
1450	end if ! of if (waterrain)
1451	end if ! of if (water)
1452
1453
1454	! 7c. Aerosol particles
1455	! -------------------
1456	! -------------
1457	! Sedimentation
1458	! -------------
1459	if (sedimentation) then
1460	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1461	zdqssed(1:ngrid,1:nq) = 0.0
1462
1463
1464	!-------------------------
1465	! find qtot
1466	if(watertest)then
1467	iq=igcm_h2o_ice
1468	call planetwide_sumval(massarea(:,:)pq(:,:,iq)ptimestep/totarea_planet,dWtot)
1469	call planetwide_sumval(massarea(:,:)pdq(:,:,iq)ptimestep/totarea_planet,dWtots)
1470	if (is_master) then
1471	print*,'Before sedim pq =',dWtot,' kg m-2'
1472	print*,'Before sedim pdq =',dWtots,' kg m-2'
1473	endif
1474	endif
1475	!-------------------------
1476
1477	call callsedim(ngrid,nlayer,ptimestep, &
1478	pplev,zzlev,pt,pdt,pq,pdq,zdqsed,zdqssed,nq)
1479
1480	!-------------------------
1481	! find qtot
1482	if(watertest)then
1483	iq=igcm_h2o_ice
1484	call planetwide_sumval(massarea(:,:)pq(:,:,iq)ptimestep/totarea_planet,dWtot)
1485	call planetwide_sumval(massarea(:,:)pdq(:,:,iq)ptimestep/totarea_planet,dWtots)
1486	if (is_master) then
1487	print*,'After sedim pq =',dWtot,' kg m-2'
1488	print*,'After sedim pdq =',dWtots,' kg m-2'
1489	endif
1490	endif
1491	!-------------------------
1492
1493	! for now, we only allow H2O ice to sediment
1494	! and as in rain.F90, whether it falls as rain or snow depends
1495	! only on the surface temperature
1496	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1497	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1498
1499	!-------------------------
1500	! test water conservation
1501	if(watertest)then
1502	call planetwide_sumval(massarea(:,:)(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))ptimestep/totarea_planet,dWtot)
1503	call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))area(:)ptimestep/totarea_planet,dWtots)
1504	if (is_master) then
1505	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
1506	print*,'In sedim surface ice change =',dWtots,' kg m-2'
1507	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
1508	endif
1509	endif
1510	!-------------------------
1511
1512	end if ! of if (sedimentation)
1513
1514
1515	! 7d. Updates
1516	! ---------
1517
1518	! -----------------------------------
1519	! Updating atm mass and tracer budget
1520	! -----------------------------------
1521
1522	if(mass_redistrib) then
1523
1524	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * &
1525	( zdqevap(1:ngrid,1:nlayer) &
1526	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) &
1527	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) &
1528	+ dqvaplscale(1:ngrid,1:nlayer) )
1529
1530	do ig = 1, ngrid
1531	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1532	enddo
1533
1534	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1535	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1536	call writediagfi(ngrid,"mass","mass"," ",3,mass)
1537
1538	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1539	rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1540	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1541	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1542
1543
1544	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1545	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1546	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1547	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1548	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1549	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1550	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1551
1552	! print,'after mass redistrib, q=',pq(211,1:nlayer,igcm_h2o_vap)+ptimesteppdq(211,1:nlayer,igcm_h2o_vap)
1553	endif
1554
1555
1556	! 7e. Ocean
1557	! ---------
1558
1559	! ---------------------------------
1560	! Slab_ocean
1561	! ---------------------------------
1562	if (ok_slab_ocean)then
1563
1564	do ig=1,ngrid
1565	qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
1566	enddo
1567
1568	call ocean_slab_ice(ptimestep, &
1569	ngrid, knindex, tsea_ice, fluxrad, &
1570	zdqssnow, qsurf(:,igcm_h2o_ice), &
1571	-zdqsdif(:,igcm_h2o_vap), &
1572	flux_sens_lat,tsea_ice, pctsrf_sic, &
1573	taux,tauy,icount)
1574
1575
1576	call ocean_slab_get_vars(ngrid,tslab, &
1577	sea_ice, flux_o, &
1578	flux_g, dt_hdiff, &
1579	dt_ekman)
1580
1581	do ig=1,ngrid
1582	if (nint(rnat(ig)).eq.1)then
1583	tslab(ig,1) = 0.
1584	tslab(ig,2) = 0.
1585	tsea_ice(ig) = 0.
1586	sea_ice(ig) = 0.
1587	pctsrf_sic(ig) = 0.
1588	qsurf(ig,igcm_h2o_ice)=qsurfint(ig,igcm_h2o_ice)
1589	endif
1590	enddo
1591
1592
1593	endif! (ok_slab_ocean)
1594
1595	! ---------------------------------
1596	! Updating tracer budget on surface
1597	! ---------------------------------
1598
1599	if(hydrology)then
1600
1601	call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
1602	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, &
1603	sea_ice)
1604	! note: for now, also changes albedo in the subroutine
1605
1606	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
1607	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqs_hyd(1:ngrid,1:nq)
1608	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
1609
1610	!-------------------------
1611	! test energy conservation
1612	if(enertest)then
1613	call planetwide_sumval(area(:)capcal(:)zdtsurf_hyd(:)/totarea_planet,dEtots)
1614	if (is_master) print*,'In hydrol surface energy change =',dEtots,' W m-2'
1615	endif
1616	!-------------------------
1617
1618	!-------------------------
1619	! test water conservation
1620	if(watertest)then
1621	call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)area(:)ptimestep/totarea_planet,dWtots)
1622	if (is_master) print*,'In hydrol surface ice change =',dWtots,' kg m-2'
1623	call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)area(:)ptimestep/totarea_planet,dWtots)
1624	if (is_master) then
1625	print*,'In hydrol surface water change =',dWtots,' kg m-2'
1626	print*,'---------------------------------------------------------------'
1627	endif
1628	endif
1629	!-------------------------
1630
1631	ELSE ! of if (hydrology)
1632
1633	qsurf(1:ngrid,1:nq)=qsurf(1:ngrid,1:nq)+ptimestep*dqsurf(1:ngrid,1:nq)
1634
1635	END IF ! of if (hydrology)
1636
1637	! Add qsurf to qsurf_hist, which is what we save in
1638	! diagfi.nc etc. At the same time, we set the water
1639	! content of ocean gridpoints back to zero, in order
1640	! to avoid rounding errors in vdifc, rain
1641	qsurf_hist(:,:) = qsurf(:,:)
1642
1643	if(ice_update)then
1644	ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
1645	endif
1646
1647	endif ! of if (tracer)
1648
1649	!-----------------------------------------------------------------------
1650	! 9. Surface and sub-surface soil temperature
1651	!-----------------------------------------------------------------------
1652
1653	! Increment surface temperature
1654
1655	if(ok_slab_ocean)then
1656	do ig=1,ngrid
1657	if (nint(rnat(ig)).eq.0)then
1658	zdtsurf(ig)= (tslab(ig,1) &
1659	+ pctsrf_sic(ig)*(tsea_ice(ig)-tslab(ig,1))-tsurf(ig))/ptimestep
1660	endif
1661	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1662	enddo
1663
1664	else
1665	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1666	endif!(ok_slab_ocean)
1667
1668	! Compute soil temperatures and subsurface heat flux
1669	if (callsoil) then
1670	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
1671	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1672	endif
1673
1674
1675	if (ok_slab_ocean) then
1676	do ig=1,ngrid
1677	fluxgrdocean(ig)=fluxgrd(ig)
1678	if (nint(rnat(ig)).eq.0) then
1679	capcal(ig)=capcalocean
1680	fluxgrd(ig)=0.
1681	fluxgrdocean(ig)=pctsrf_sic(ig)flux_g(ig)+(1-pctsrf_sic(ig))(dt_hdiff(ig,1)+dt_ekman(ig,1))
1682	do iq=1,nsoilmx
1683	tsoil(ig,iq)=tsurf(ig)
1684	enddo
1685	if (pctsrf_sic(ig).gt.0.5) then
1686	capcal(ig)=capcalseaice
1687	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
1688	capcal(ig)=capcalsno
1689	endif
1690	endif
1691	endif
1692	enddo
1693	endif !(ok_slab_ocean)
1694
1695	!-------------------------
1696	! test energy conservation
1697	if(enertest)then
1698	call planetwide_sumval(area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1699	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1700	endif
1701	!-------------------------
1702
1703	!-----------------------------------------------------------------------
1704	! 10. Perform diagnostics and write output files
1705	!-----------------------------------------------------------------------
1706
1707	! -------------------------------
1708	! Dynamical fields incrementation
1709	! -------------------------------
1710	! For output only: the actual model integration is performed in the dynamics
1711
1712	! temperature, zonal and meridional wind
1713	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1714	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1715	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1716
1717	! diagnostic
1718	zdtdyn(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)
1719	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1720
1721	if(firstcall)then
1722	zdtdyn(1:ngrid,1:nlayer)=0.0
1723	endif
1724
1725	! dynamical heating diagnostic
1726	do ig=1,ngrid
1727	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp/ptimestep
1728	enddo
1729
1730	! tracers
1731	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1732
1733	! surface pressure
1734	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1735
1736	! pressure
1737	do l=1,nlayer
1738	zplev(1:ngrid,l) = pplev(1:ngrid,l)/pplev(1:ngrid,1)*ps(:)
1739	zplay(1:ngrid,l) = pplay(1:ngrid,l)/pplev(1:ngrid,1)*ps(1:ngrid)
1740	enddo
1741
1742	! ---------------------------------------------------------
1743	! Surface and soil temperature information
1744	! ---------------------------------------------------------
1745
1746	call planetwide_sumval(area(:)*tsurf(:)/totarea_planet,Ts1)
1747	call planetwide_minval(tsurf(:),Ts2)
1748	call planetwide_maxval(tsurf(:),Ts3)
1749	if(callsoil)then
1750	TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1751	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1752	print*,Ts1,Ts2,Ts3,TsS
1753	else
1754	if (is_master) then
1755	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1756	print*,Ts1,Ts2,Ts3
1757	endif
1758	end if
1759
1760	! ---------------------------------------------------------
1761	! Check the energy balance of the simulation during the run
1762	! ---------------------------------------------------------
1763
1764	if(corrk)then
1765
1766	call planetwide_sumval(area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1767	call planetwide_sumval(area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1768	call planetwide_sumval(area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1769	call planetwide_sumval(area(:)*fluxgrd(:)/totarea_planet,GND)
1770	call planetwide_sumval(area(:)*fluxdyn(:)/totarea_planet,DYN)
1771	do ig=1,ngrid
1772	if(fluxtop_dn(ig).lt.0.0)then
1773	print*,'fluxtop_dn has gone crazy'
1774	print*,'fluxtop_dn=',fluxtop_dn(ig)
1775	print*,'tau_col=',tau_col(ig)
1776	print*,'aerosol=',aerosol(ig,:,:)
1777	print*,'temp= ',pt(ig,:)
1778	print*,'pplay= ',pplay(ig,:)
1779	call abort
1780	endif
1781	end do
1782
1783	if(ngrid.eq.1)then
1784	DYN=0.0
1785	endif
1786
1787	if (is_master) then
1788	print*,' ISR ASR OLR GND DYN [W m^-2]'
1789	print*, ISR,ASR,OLR,GND,DYN
1790	endif
1791
1792	if(enertest .and. is_master)then
1793	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1794	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1795	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1796	endif
1797
1798	if(meanOLR .and. is_master)then
1799	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1800	! to record global radiative balance
1801	open(92,file="rad_bal.out",form='formatted',position='append')
1802	write(92,*) zday,ISR,ASR,OLR
1803	close(92)
1804	open(93,file="tem_bal.out",form='formatted',position='append')
1805	if(callsoil)then
1806	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1807	else
1808	write(93,*) zday,Ts1,Ts2,Ts3
1809	endif
1810	close(93)
1811	endif
1812	endif
1813
1814	endif
1815
1816
1817	! ------------------------------------------------------------------
1818	! Diagnostic to test radiative-convective timescales in code
1819	! ------------------------------------------------------------------
1820	if(testradtimes)then
1821	open(38,file="tau_phys.out",form='formatted',position='append')
1822	ig=1
1823	do l=1,nlayer
1824	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1825	enddo
1826	close(38)
1827	print*,'As testradtimes enabled,'
1828	print*,'exiting physics on first call'
1829	call abort
1830	endif
1831
1832	! ---------------------------------------------------------
1833	! Compute column amounts (kg m-2) if tracers are enabled
1834	! ---------------------------------------------------------
1835	if(tracer)then
1836	qcol(1:ngrid,1:nq)=0.0
1837	do iq=1,nq
1838	do ig=1,ngrid
1839	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1840	enddo
1841	enddo
1842
1843	! Generalised for arbitrary aerosols now. (LK)
1844	reffcol(1:ngrid,1:naerkind)=0.0
1845	if(co2cond.and.(iaero_co2.ne.0))then
1846	call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2))
1847	do ig=1,ngrid
1848	reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)reffrad(ig,1:nlayer,iaero_co2)mass(ig,1:nlayer))
1849	enddo
1850	endif
1851	if(water.and.(iaero_h2o.ne.0))then
1852	call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, &
1853	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
1854	do ig=1,ngrid
1855	reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)reffrad(ig,1:nlayer,iaero_h2o)mass(ig,1:nlayer))
1856	enddo
1857	endif
1858
1859	endif
1860
1861	! ---------------------------------------------------------
1862	! Test for water conservation if water is enabled
1863	! ---------------------------------------------------------
1864
1865	if(water)then
1866
1867	call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf)
1868	call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf)
1869	call planetwide_sumval(area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol)
1870	call planetwide_sumval(area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol)
1871
1872	h2otot = icesrf + liqsrf + icecol + vapcol
1873
1874	if (is_master) then
1875	print*,' Total water amount [kg m^-2]: ',h2otot
1876	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
1877	print*, icesrf,liqsrf,icecol,vapcol
1878	endif
1879
1880	if(meanOLR .and. is_master)then
1881	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1882	! to record global water balance
1883	open(98,file="h2o_bal.out",form='formatted',position='append')
1884	write(98,*) zday,icesrf,liqsrf,icecol,vapcol
1885	close(98)
1886	endif
1887	endif
1888
1889	endif
1890
1891	! ---------------------------------------------------------
1892	! Calculate RH for diagnostic if water is enabled
1893	! ---------------------------------------------------------
1894
1895	if(water)then
1896	do l = 1, nlayer
1897	do ig=1,ngrid
1898	call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
1899	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
1900	enddo
1901	enddo
1902
1903	! compute maximum possible H2O column amount (100% saturation)
1904	do ig=1,ngrid
1905	H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
1906	enddo
1907
1908	endif
1909
1910
1911	print,'--> Ls =',zls180./pi
1912	! -------------------------------------------------------------------
1913	! Writing NetCDF file "RESTARTFI" at the end of the run
1914	! -------------------------------------------------------------------
1915	! Note: 'restartfi' is stored just before dynamics are stored
1916	! in 'restart'. Between now and the writting of 'restart',
1917	! there will have been the itau=itau+1 instruction and
1918	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1919	! thus we store for time=time+dtvr
1920
1921	if(lastcall) then
1922	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1923
1924
1925	! Update surface ice distribution to iterate to steady state if requested
1926	if(ice_update)then
1927
1928	do ig=1,ngrid
1929
1930	delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
1931
1932	! add multiple years of evolution
1933	qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
1934
1935	! if ice has gone -ve, set to zero
1936	if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
1937	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1938	endif
1939
1940	! if ice is seasonal, set to zero (NEW)
1941	if(ice_min(ig).lt.0.01)then
1942	qsurf_hist(ig,igcm_h2o_ice) = 0.0
1943	endif
1944
1945	enddo
1946
1947	! enforce ice conservation
1948	ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
1949	qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
1950
1951	endif
1952
1953	if (ngrid.ne.1) then
1954	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1955	!#ifdef CPP_PARA
1956	!! for now in parallel we use a different routine to write restartfi.nc
1957	! call phyredem(ngrid,"restartfi.nc", &
1958	! ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1959	! cloudfrac,totcloudfrac,hice)
1960	!#else
1961	! call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq, &
1962	! ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1963	! area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1964	! cloudfrac,totcloudfrac,hice,noms)
1965	!#endif
1966	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
1967	ptimestep,ztime_fin, &
1968	tsurf,tsoil,emis,q2,qsurf_hist, &
1969	cloudfrac,totcloudfrac,hice, &
1970	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1971	endif
1972
1973	if(ok_slab_ocean) then
1974	call ocean_slab_final!(tslab, seaice)
1975	end if
1976
1977
1978	endif ! of if (lastcall)
1979
1980
1981	! -----------------------------------------------------------------
1982	! Saving statistics :
1983	! -----------------------------------------------------------------
1984	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1985	! which can later be used to make the statistic files of the run:
1986	! "stats") only possible in 3D runs !
1987
1988
1989	if (callstats) then
1990
1991	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1992	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1993	call wstats(ngrid,"fluxsurf_lw", &
1994	"Thermal IR radiative flux to surface","W.m-2",2, &
1995	fluxsurf_lw)
1996	! call wstats(ngrid,"fluxsurf_sw", &
1997	! "Solar radiative flux to surface","W.m-2",2, &
1998	! fluxsurf_sw_tot)
1999	call wstats(ngrid,"fluxtop_lw", &
2000	"Thermal IR radiative flux to space","W.m-2",2, &
2001	fluxtop_lw)
2002	! call wstats(ngrid,"fluxtop_sw", &
2003	! "Solar radiative flux to space","W.m-2",2, &
2004	! fluxtop_sw_tot)
2005
2006	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2007	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2008	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2009	call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
2010	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
2011	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2012	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2013	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
2014	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
2015	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
2016
2017	if (tracer) then
2018	do iq=1,nq
2019	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2020	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2021	'kg m^-2',2,qsurf(1,iq) )
2022
2023	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2024	'kg m^-2',2,qcol(1,iq) )
2025	! call wstats(ngrid,trim(noms(iq))//'_reff', &
2026	! trim(noms(iq))//'_reff', &
2027	! 'm',3,reffrad(1,1,iq))
2028	end do
2029	if (water) then
2030	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
2031	call wstats(ngrid,"vmr_h2ovapor", &
2032	"H2O vapour volume mixing ratio","mol/mol", &
2033	3,vmr)
2034	endif ! of if (water)
2035
2036	endif !tracer
2037
2038	if(watercond.and.CLFvarying)then
2039	call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
2040	call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
2041	call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2042	call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2043	call wstats(ngrid,"RH","relative humidity"," ",3,RH)
2044	endif
2045
2046
2047
2048	if (ok_slab_ocean) then
2049	call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
2050	call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
2051	call wstats(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
2052	call wstats(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
2053	call wstats(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
2054	call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
2055	call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
2056	call wstats(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
2057	call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
2058	call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
2059	endif! (ok_slab_ocean)
2060
2061	if(lastcall) then
2062	write (,) "Writing stats..."
2063	call mkstats(ierr)
2064	endif
2065	endif !if callstats
2066
2067
2068	! ----------------------------------------------------------------------
2069	! output in netcdf file "DIAGFI", containing any variable for diagnostic
2070	! (output with period "ecritphy", set in "run.def")
2071	! ----------------------------------------------------------------------
2072	! writediagfi can also be called from any other subroutine for any variable.
2073	! but its preferable to keep all the calls in one place...
2074
2075	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
2076	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2077	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2078	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2079	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
2080	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2081	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2082	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2083	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
2084
2085	! Subsurface temperatures
2086	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2087	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
2088
2089	! Oceanic layers
2090	if(ok_slab_ocean) then
2091	call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
2092	call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
2093	call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
2094	call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
2095	call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
2096	call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
2097	call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
2098	call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
2099	call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
2100	call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
2101	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2102	call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
2103	endif! (ok_slab_ocean)
2104
2105	! Total energy balance diagnostics
2106	if(callrad.and.(.not.newtonian))then
2107	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
2108	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2109	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2110	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2111	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
2112	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
2113	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
2114	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
2115	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
2116	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
2117	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
2118	if(ok_slab_ocean) then
2119	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
2120	else
2121	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2122	endif!(ok_slab_ocean)
2123	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2124	endif
2125
2126	if(enertest) then
2127	if (calldifv) then
2128	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
2129	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
2130	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
2131	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
2132	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2133	endif
2134	if (corrk) then
2135	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
2136	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
2137	endif
2138	if(watercond) then
2139	! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
2140	! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
2141	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
2142	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
2143	call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
2144	! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
2145	endif
2146	endif
2147
2148	! Temporary inclusions for heating diagnostics
2149	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2150	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2151	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2152	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2153
2154	! debugging
2155	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2156	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2157
2158	! Output aerosols
2159	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
2160	call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
2161	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
2162	call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
2163	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
2164	call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
2165	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
2166	call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
2167
2168	! Output tracers
2169	if (tracer) then
2170	do iq=1,nq
2171	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2172	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2173	! 'kg m^-2',2,qsurf(1,iq) )
2174	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2175	'kg m^-2',2,qsurf_hist(1,iq) )
2176	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2177	'kg m^-2',2,qcol(1,iq) )
2178
2179	if(watercond.or.CLFvarying)then
2180	call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
2181	call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
2182	call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2183	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2184	call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
2185	endif
2186
2187	if(waterrain)then
2188	call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
2189	call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
2190	endif
2191
2192	if((hydrology).and.(.not.ok_slab_ocean))then
2193	call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
2194	endif
2195
2196	if(ice_update)then
2197	call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
2198	call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
2199	endif
2200
2201	do ig=1,ngrid
2202	if(tau_col(ig).gt.1.e3)then
2203	! print*,'WARNING: tau_col=',tau_col(ig)
2204	! print*,'at ig=',ig,'in PHYSIQ'
2205	endif
2206	end do
2207
2208	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2209
2210	enddo
2211	endif
2212
2213	! output spectrum
2214	if(specOLR.and.corrk)then
2215	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2216	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2217	endif
2218
2219	icount=icount+1
2220
2221	if (lastcall) then
2222
2223	! deallocate gas variables
2224	!$OMP BARRIER
2225	!$OMP MASTER
2226	IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
2227	IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac ) ! both allocated in su_gases.F90
2228	!$OMP END MASTER
2229	!$OMP BARRIER
2230
2231	! deallocate saved arrays
2232	! this is probably not that necessary
2233	! ... but probably a good idea to clean the place before we leave
2234	IF ( ALLOCATED(tsurf)) DEALLOCATE(tsurf)
2235	IF ( ALLOCATED(tsoil)) DEALLOCATE(tsoil)
2236	IF ( ALLOCATED(albedo)) DEALLOCATE(albedo)
2237	IF ( ALLOCATED(albedo0)) DEALLOCATE(albedo0)
2238	IF ( ALLOCATED(rnat)) DEALLOCATE(rnat)
2239	IF ( ALLOCATED(emis)) DEALLOCATE(emis)
2240	IF ( ALLOCATED(dtrad)) DEALLOCATE(dtrad)
2241	IF ( ALLOCATED(fluxrad_sky)) DEALLOCATE(fluxrad_sky)
2242	IF ( ALLOCATED(fluxrad)) DEALLOCATE(fluxrad)
2243	IF ( ALLOCATED(capcal)) DEALLOCATE(capcal)
2244	IF ( ALLOCATED(fluxgrd)) DEALLOCATE(fluxgrd)
2245	IF ( ALLOCATED(qsurf)) DEALLOCATE(qsurf)
2246	IF ( ALLOCATED(q2)) DEALLOCATE(q2)
2247	IF ( ALLOCATED(ztprevious)) DEALLOCATE(ztprevious)
2248	IF ( ALLOCATED(hice)) DEALLOCATE(hice)
2249	IF ( ALLOCATED(cloudfrac)) DEALLOCATE(cloudfrac)
2250	IF ( ALLOCATED(totcloudfrac)) DEALLOCATE(totcloudfrac)
2251	IF ( ALLOCATED(qsurf_hist)) DEALLOCATE(qsurf_hist)
2252	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
2253	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
2254	IF ( ALLOCATED(ice_initial)) DEALLOCATE(ice_initial)
2255	IF ( ALLOCATED(ice_min)) DEALLOCATE(ice_min)
2256
2257	IF ( ALLOCATED(fluxsurf_lw)) DEALLOCATE(fluxsurf_lw)
2258	IF ( ALLOCATED(fluxsurf_sw)) DEALLOCATE(fluxsurf_sw)
2259	IF ( ALLOCATED(fluxtop_lw)) DEALLOCATE(fluxtop_lw)
2260	IF ( ALLOCATED(fluxabs_sw)) DEALLOCATE(fluxabs_sw)
2261	IF ( ALLOCATED(fluxtop_dn)) DEALLOCATE(fluxtop_dn)
2262	IF ( ALLOCATED(fluxdyn)) DEALLOCATE(fluxdyn)
2263	IF ( ALLOCATED(OLR_nu)) DEALLOCATE(OLR_nu)
2264	IF ( ALLOCATED(OSR_nu)) DEALLOCATE(OSR_nu)
2265	IF ( ALLOCATED(sensibFlux)) DEALLOCATE(sensibFlux)
2266	IF ( ALLOCATED(zdtlw)) DEALLOCATE(zdtlw)
2267	IF ( ALLOCATED(zdtsw)) DEALLOCATE(zdtsw)
2268	IF ( ALLOCATED(tau_col)) DEALLOCATE(tau_col)
2269	IF ( ALLOCATED(pctsrf_sic)) DEALLOCATE(pctsrf_sic)
2270	IF ( ALLOCATED(tslab)) DEALLOCATE(tslab)
2271	IF ( ALLOCATED(tsea_ice)) DEALLOCATE(tsea_ice)
2272	IF ( ALLOCATED(sea_ice)) DEALLOCATE(sea_ice)
2273	IF ( ALLOCATED(zmasq)) DEALLOCATE(zmasq)
2274	IF ( ALLOCATED(knindex)) DEALLOCATE(knindex)
2275
2276	!! this is defined in comsaison_h
2277	IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
2278
2279	IF ( ALLOCATED(fract)) DEALLOCATE(fract)
2280
2281
2282	!! this is defined in radcommon_h
2283	IF ( ALLOCATED(eclipse)) DEALLOCATE(eclipse)
2284
2285	!! this is defined in comsoil_h
2286	IF ( ALLOCATED(layer)) DEALLOCATE(layer)
2287	IF ( ALLOCATED(mlayer)) DEALLOCATE(mlayer)
2288	IF ( ALLOCATED(inertiedat)) DEALLOCATE(inertiedat)
2289
2290	!! this is defined in surfdat_h
2291	IF ( ALLOCATED(albedodat)) DEALLOCATE(albedodat)
2292	IF ( ALLOCATED(phisfi)) DEALLOCATE(phisfi)
2293	IF ( ALLOCATED(zmea)) DEALLOCATE(zmea)
2294	IF ( ALLOCATED(zstd)) DEALLOCATE(zstd)
2295	IF ( ALLOCATED(zsig)) DEALLOCATE(zsig)
2296	IF ( ALLOCATED(zgam)) DEALLOCATE(zgam)
2297	IF ( ALLOCATED(zthe)) DEALLOCATE(zthe)
2298	IF ( ALLOCATED(dryness)) DEALLOCATE(dryness)
2299	IF ( ALLOCATED(watercaptag)) DEALLOCATE(watercaptag)
2300
2301	!! this is defined in tracer_h
2302	IF ( ALLOCATED(noms)) DEALLOCATE(noms)
2303	IF ( ALLOCATED(mmol)) DEALLOCATE(mmol)
2304	IF ( ALLOCATED(radius)) DEALLOCATE(radius)
2305	IF ( ALLOCATED(rho_q)) DEALLOCATE(rho_q)
2306	IF ( ALLOCATED(qext)) DEALLOCATE(qext)
2307	IF ( ALLOCATED(alpha_lift)) DEALLOCATE(alpha_lift)
2308	IF ( ALLOCATED(alpha_devil)) DEALLOCATE(alpha_devil)
2309	IF ( ALLOCATED(qextrhor)) DEALLOCATE(qextrhor)
2310	IF ( ALLOCATED(igcm_dustbin)) DEALLOCATE(igcm_dustbin)
2311
2312	!! this is defined in comgeomfi_h
2313	IF ( ALLOCATED(lati)) DEALLOCATE(lati)
2314	IF ( ALLOCATED(long)) DEALLOCATE(long)
2315	IF ( ALLOCATED(area)) DEALLOCATE(area)
2316	IF ( ALLOCATED(sinlat)) DEALLOCATE(sinlat)
2317	IF ( ALLOCATED(coslat)) DEALLOCATE(coslat)
2318	IF ( ALLOCATED(sinlon)) DEALLOCATE(sinlon)
2319	IF ( ALLOCATED(coslon)) DEALLOCATE(coslon)
2320	endif
2321
2322	return
2323	end subroutine physiq

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

Download in other formats: