Context Navigation

physiq.F90 @ 1371

Last change on this file since 1371 was 1346, checked in by emillour, 11 years ago

Mars and Generic GCM:

Bug correction on the computation of vertical velocity.

Property svn:executable set to *

File size: 94.4 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, right_ascen
20	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
21	USE comgeomfi_h, only: long, lati, area, totarea, totarea_planet
22	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
23	alpha_lift, alpha_devil, qextrhor, &
24	igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, &
25	igcm_co2_ice
26	use control_mod, only: ecritphy, iphysiq, nday
27	use phyredem, only: physdem0, physdem1
28	use slab_ice_h
29	use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
30	ocean_slab_get_vars,ocean_slab_final
31	use surf_heat_transp_mod,only :init_masquv
32	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
33	use mod_phys_lmdz_para, only : is_master
34	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
35	obliquit, nres, z0
36
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 zlss ! sub solar point longitude (rad)
238	real zday ! date (time since Ls=0, in martian days)
239	real zzlay(ngrid,nlayer) ! altitude at the middle of the layers
240	real zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
241	real latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
242
243	! Tendencies due to various processes:
244	real dqsurf(ngrid,nq)
245	real cldtlw(ngrid,nlayer) ! (K/s) LW heating rate for clear areas
246	real cldtsw(ngrid,nlayer) ! (K/s) SW heating rate for clear areas
247	real zdtsurf(ngrid) ! (K/s)
248	real dtlscale(ngrid,nlayer)
249	real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
250	real zdhdif(ngrid,nlayer), zdtsdif(ngrid) ! (K/s)
251	real zdtdif(ngrid,nlayer) ! (K/s)
252	real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
253	real zdhadj(ngrid,nlayer) ! (K/s)
254	real zdtgw(ngrid,nlayer) ! (K/s)
255	real zdtmr(ngrid,nlayer)
256	real zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
257	real zdtc(ngrid,nlayer),zdtsurfc(ngrid)
258	real zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
259	real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)
260	real zdtsurfmr(ngrid)
261
262	real zdmassmr(ngrid,nlayer),zdpsrfmr(ngrid)
263	real zdmassmr_col(ngrid)
264
265	real zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
266	real zdqevap(ngrid,nlayer)
267	real zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
268	real zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
269	real zdqadj(ngrid,nlayer,nq)
270	real zdqc(ngrid,nlayer,nq)
271	real zdqmr(ngrid,nlayer,nq),zdqsurfmr(ngrid,nq)
272	real zdqlscale(ngrid,nlayer,nq)
273	real zdqslscale(ngrid,nq)
274	real zdqchim(ngrid,nlayer,nq)
275	real zdqschim(ngrid,nq)
276
277	real zdteuv(ngrid,nlayer) ! (K/s)
278	real zdtconduc(ngrid,nlayer) ! (K/s)
279	real zdumolvis(ngrid,nlayer)
280	real zdvmolvis(ngrid,nlayer)
281	real zdqmoldiff(ngrid,nlayer,nq)
282
283	! Local variables for local calculations:
284	real zflubid(ngrid)
285	real zplanck(ngrid),zpopsk(ngrid,nlayer)
286	real zdum1(ngrid,nlayer)
287	real zdum2(ngrid,nlayer)
288	real ztim1,ztim2,ztim3, z1,z2
289	real ztime_fin
290	real zdh(ngrid,nlayer)
291	real gmplanet
292	real taux(ngrid),tauy(ngrid)
293
294	integer length
295	parameter (length=100)
296
297	! local variables only used for diagnostics (output in file "diagfi" or "stats")
298	! ------------------------------------------------------------------------------
299	real ps(ngrid), zt(ngrid,nlayer)
300	real zu(ngrid,nlayer),zv(ngrid,nlayer)
301	real zq(ngrid,nlayer,nq)
302	character*2 str2
303	character*5 str5
304	real zdtadj(ngrid,nlayer)
305	real zdtdyn(ngrid,nlayer)
306	real,allocatable,dimension(:,:),save :: ztprevious
307	!$OMP THREADPRIVATE(ztprevious)
308	real reff(ngrid,nlayer) ! effective dust radius (used if doubleq=T)
309	real qtot1,qtot2 ! total aerosol mass
310	integer igmin, lmin
311	logical tdiag
312
313	real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
314	real zstress(ngrid), cd
315	real hco2(nq), tmean, zlocal(nlayer)
316	real vmr(ngrid,nlayer) ! volume mixing ratio
317
318	real time_phys
319
320	! reinstated by RW for diagnostic
321	real,allocatable,dimension(:),save :: tau_col
322	!$OMP THREADPRIVATE(tau_col)
323
324	! included by RW to reduce insanity of code
325	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
326
327	! included by RW to compute tracer column densities
328	real qcol(ngrid,nq)
329
330	! included by RW for H2O precipitation
331	real zdtrain(ngrid,nlayer)
332	real zdqrain(ngrid,nlayer,nq)
333	real zdqsrain(ngrid)
334	real zdqssnow(ngrid)
335	real rainout(ngrid)
336
337	! included by RW for H2O Manabe scheme
338	real dtmoist(ngrid,nlayer)
339	real dqmoist(ngrid,nlayer,nq)
340
341	real qvap(ngrid,nlayer)
342	real dqvaplscale(ngrid,nlayer)
343	real dqcldlscale(ngrid,nlayer)
344	real rneb_man(ngrid,nlayer)
345	real rneb_lsc(ngrid,nlayer)
346
347	! included by RW to account for surface cooling by evaporation
348	real dtsurfh2olat(ngrid)
349
350
351	! to test energy conservation (RW+JL)
352	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
353	real dEtot, dEtots, AtmToSurf_TurbFlux
354	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
355	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
356	real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
357	real dEdiffs(ngrid),dEdiff(ngrid)
358	real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayer), lscaledEz(ngrid,nlayer)
359	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
360	real dtmoist_max,dtmoist_min
361
362	real dItot, dItot_tmp, dVtot, dVtot_tmp
363
364	! included by BC for evaporation
365	real qevap(ngrid,nlayer,nq)
366	real tevap(ngrid,nlayer)
367	real dqevap1(ngrid,nlayer)
368	real dtevap1(ngrid,nlayer)
369
370	! included by BC for hydrology
371	real,allocatable,save :: hice(:)
372	!$OMP THREADPRIVATE(hice)
373
374	! included by RW to test water conservation (by routine)
375	real h2otot
376	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
377	real h2o_surf_all
378	logical watertest
379	save watertest
380	!$OMP THREADPRIVATE(watertest)
381
382	! included by RW for RH diagnostic
383	real qsat(ngrid,nlayer), RH(ngrid,nlayer), H2Omaxcol(ngrid),psat_tmp
384
385	! included by RW for hydrology
386	real dqs_hyd(ngrid,nq)
387	real zdtsurf_hyd(ngrid)
388
389	! included by RW for water cycle conservation tests
390	real icesrf,liqsrf,icecol,vapcol
391
392	! included by BC for double radiative transfer call
393	logical clearsky
394	real zdtsw1(ngrid,nlayer)
395	real zdtlw1(ngrid,nlayer)
396	real fluxsurf_lw1(ngrid)
397	real fluxsurf_sw1(ngrid)
398	real fluxtop_lw1(ngrid)
399	real fluxabs_sw1(ngrid)
400	real tau_col1(ngrid)
401	real OLR_nu1(ngrid,L_NSPECTI)
402	real OSR_nu1(ngrid,L_NSPECTV)
403	real tf, ntf
404
405	! included by BC for cloud fraction computations
406	real,allocatable,dimension(:,:),save :: cloudfrac
407	real,allocatable,dimension(:),save :: totcloudfrac
408	!$OMP THREADPRIVATE(cloudfrac,totcloudfrac)
409
410	! included by RW for vdifc water conservation test
411	real nconsMAX
412	real vdifcncons(ngrid)
413	real cadjncons(ngrid)
414
415	! double precision qsurf_hist(ngrid,nq)
416	real,allocatable,dimension(:,:),save :: qsurf_hist
417	!$OMP THREADPRIVATE(qsurf_hist)
418
419	! included by RW for temp convadj conservation test
420	real playtest(nlayer)
421	real plevtest(nlayer)
422	real ttest(nlayer)
423	real qtest(nlayer)
424	integer igtest
425
426	! included by RW for runway greenhouse 1D study
427	real muvar(ngrid,nlayer+1)
428
429	! included by RW for variable H2O particle sizes
430	real,allocatable,dimension(:,:,:),save :: reffrad ! aerosol effective radius (m)
431	real,allocatable,dimension(:,:,:),save :: nueffrad ! aerosol effective radius variance
432	!$OMP THREADPRIVATE(reffrad,nueffrad)
433	! real :: nueffrad_dummy(ngrid,nlayer,naerkind) !! AS. This is temporary. Check below why.
434	real :: reffh2oliq(ngrid,nlayer) ! liquid water particles effective radius (m)
435	real :: reffh2oice(ngrid,nlayer) ! water ice particles effective radius (m)
436	! real reffH2O(ngrid,nlayer)
437	real reffcol(ngrid,naerkind)
438
439	! included by RW for sourceevol
440	real,allocatable,dimension(:),save :: ice_initial
441	real delta_ice,ice_tot
442	real,allocatable,dimension(:),save :: ice_min
443
444	integer num_run
445	logical,save :: ice_update
446	!$OMP THREADPRIVATE(ice_initial,ice_min,ice_update)
447
448	! included by MS to compute the daily average of rings shadowing
449	integer, parameter :: nb_hours = 1536 ! set how many times per day are used
450	real :: pas
451	integer :: m
452	! temporary variables computed at each step of this average
453	real :: ptime_day ! Universal time in sol fraction
454	real :: tmp_zls,tmp_dist_star, tmp_declin, tmp_right_ascen ! tmp solar longitude, stellar dist, declin and RA
455	real, dimension(:), allocatable :: tmp_fract ! day fraction of the time interval
456	real, dimension(:), allocatable :: tmp_mu0 ! equivalent solar angle
457
458	! included by BC for slab ocean
459	real, dimension(:),allocatable,save :: pctsrf_sic
460	real, dimension(:,:),allocatable,save :: tslab
461	real, dimension(:),allocatable,save :: tsea_ice
462	real, dimension(:),allocatable,save :: sea_ice
463	real, dimension(:),allocatable,save :: zmasq
464	integer, dimension(:),allocatable,save ::knindex
465	!$OMP THREADPRIVATE(pctsrf_sic,tslab,tsea_ice,sea_ice,zmasq,knindex)
466
467	real :: tsurf2(ngrid)
468	real :: flux_o(ngrid),flux_g(ngrid),fluxgrdocean(ngrid)
469	real :: dt_ekman(ngrid,noceanmx),dt_hdiff(ngrid,noceanmx)
470	real :: flux_sens_lat(ngrid)
471	real :: qsurfint(ngrid,nq)
472
473
474	!=======================================================================
475
476	! 1. Initialisation
477	! -----------------
478
479	! 1.1 Initialisation only at first call
480	! ---------------------------------------
481	if (firstcall) then
482
483	!!! ALLOCATE SAVED ARRAYS
484	ALLOCATE(tsurf(ngrid))
485	ALLOCATE(tsoil(ngrid,nsoilmx))
486	ALLOCATE(albedo(ngrid))
487	ALLOCATE(albedo0(ngrid))
488	ALLOCATE(rnat(ngrid))
489	ALLOCATE(emis(ngrid))
490	ALLOCATE(dtrad(ngrid,nlayer))
491	ALLOCATE(fluxrad_sky(ngrid))
492	ALLOCATE(fluxrad(ngrid))
493	ALLOCATE(capcal(ngrid))
494	ALLOCATE(fluxgrd(ngrid))
495	ALLOCATE(qsurf(ngrid,nq))
496	ALLOCATE(q2(ngrid,nlayer+1))
497	ALLOCATE(ztprevious(ngrid,nlayer))
498	ALLOCATE(cloudfrac(ngrid,nlayer))
499	ALLOCATE(totcloudfrac(ngrid))
500	ALLOCATE(hice(ngrid))
501	ALLOCATE(qsurf_hist(ngrid,nq))
502	ALLOCATE(reffrad(ngrid,nlayer,naerkind))
503	ALLOCATE(nueffrad(ngrid,nlayer,naerkind))
504	ALLOCATE(ice_initial(ngrid))
505	ALLOCATE(ice_min(ngrid))
506	ALLOCATE(fluxsurf_lw(ngrid))
507	ALLOCATE(fluxsurf_sw(ngrid))
508	ALLOCATE(fluxtop_lw(ngrid))
509	ALLOCATE(fluxabs_sw(ngrid))
510	ALLOCATE(fluxtop_dn(ngrid))
511	ALLOCATE(fluxdyn(ngrid))
512	ALLOCATE(OLR_nu(ngrid,L_NSPECTI))
513	ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
514	ALLOCATE(sensibFlux(ngrid))
515	ALLOCATE(zdtlw(ngrid,nlayer))
516	ALLOCATE(zdtsw(ngrid,nlayer))
517	ALLOCATE(tau_col(ngrid))
518	ALLOCATE(pctsrf_sic(ngrid))
519	ALLOCATE(tslab(ngrid,noceanmx))
520	ALLOCATE(tsea_ice(ngrid))
521	ALLOCATE(sea_ice(ngrid))
522	ALLOCATE(zmasq(ngrid))
523	ALLOCATE(knindex(ngrid))
524	! ALLOCATE(qsurfint(ngrid,nqmx))
525
526
527	!! this is defined in comsaison_h
528	ALLOCATE(mu0(ngrid))
529	ALLOCATE(fract(ngrid))
530
531
532	!! this is defined in radcommon_h
533	ALLOCATE(eclipse(ngrid))
534	ALLOCATE(glat(ngrid))
535
536	! variables set to 0
537	! ~~~~~~~~~~~~~~~~~~
538	dtrad(:,:) = 0.0
539	fluxrad(:) = 0.0
540	tau_col(:) = 0.0
541	zdtsw(:,:) = 0.0
542	zdtlw(:,:) = 0.0
543
544	! initialize aerosol indexes
545	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
546	call iniaerosol()
547
548
549	! initialize tracer names, indexes and properties
550	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
551	tracerdyn=tracer
552	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) !! because noms is an argument of physdem1
553	!! whether or not tracer is on
554	if (tracer) then
555	call initracer(ngrid,nq,nametrac)
556	endif ! end tracer
557
558	!
559
560	! read startfi (initial parameters)
561	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
562	call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
563	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
564	cloudfrac,totcloudfrac,hice, &
565	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
566
567	if (pday.ne.day_ini) then
568	write(,) "ERROR in physiq.F90:"
569	write(,) "bad synchronization between physics and dynamics"
570	write(,) "dynamics day: ",pday
571	write(,) "physics day: ",day_ini
572	stop
573	endif
574
575	write (,) 'In physiq day_ini =', day_ini
576
577	! Initialize albedo and orbital calculation
578	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
579	call surfini(ngrid,nq,qsurf,albedo0)
580	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
581
582	albedo(:)=albedo0(:)
583
584	if(tlocked)then
585	print*,'Planet is tidally locked at resonance n=',nres
586	print*,'Make sure you have the right rotation rate!!!'
587	endif
588
589	! Initialize oceanic variables
590	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
591
592	if (ok_slab_ocean)then
593
594	call ocean_slab_init(ngrid,ptimestep, tslab, &
595	sea_ice, pctsrf_sic)
596
597	knindex(:) = 0
598
599	do ig=1,ngrid
600	zmasq(ig)=1
601	knindex(ig) = ig
602	if (nint(rnat(ig)).eq.0) then
603	zmasq(ig)=0
604	endif
605	enddo
606
607
608	CALL init_masquv(ngrid,zmasq)
609
610
611	endif
612
613
614	! initialize soil
615	! ~~~~~~~~~~~~~~~
616	if (callsoil) then
617	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
618	ptimestep,tsurf,tsoil,capcal,fluxgrd)
619
620	if (ok_slab_ocean) then
621	do ig=1,ngrid
622	if (nint(rnat(ig)).eq.2) then
623	capcal(ig)=capcalocean
624	if (pctsrf_sic(ig).gt.0.5) then
625	capcal(ig)=capcalseaice
626	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
627	capcal(ig)=capcalsno
628	endif
629	endif
630	endif
631	enddo
632	endif
633
634
635
636
637	else
638	print*,'WARNING! Thermal conduction in the soil turned off'
639	capcal(:)=1.e6
640	fluxgrd(:)=intheat
641	print*,'Flux from ground = ',intheat,' W m^-2'
642	endif
643	icount=1
644
645	! decide whether to update ice at end of run
646	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
647	ice_update=.false.
648	if(sourceevol)then
649	!$OMP MASTER
650	open(128,file='num_run',form='formatted', &
651	status="old",iostat=ierr)
652	if (ierr.ne.0) then
653	write(,) "physiq: Error! No num_run file!"
654	write(,) " (which is needed for sourceevol option)"
655	stop
656	endif
657	read(128,*) num_run
658	close(128)
659	!$OMP END MASTER
660	!$OMP BARRIER
661
662	if(num_run.ne.0.and.mod(num_run,2).eq.0)then
663	!if(num_run.ne.0.and.mod(num_run,3).eq.0)then
664	print*,'Updating ice at end of this year!'
665	ice_update=.true.
666	ice_min(:)=1.e4
667	endif
668	endif
669
670	! In newstart now, will have to be remove (BC 2014)
671	! define surface as continent or ocean
672	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
673	if (.not.ok_slab_ocean) then
674	rnat(:)=1.
675	do ig=1,ngrid
676	! if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
677	if(inertiedat(ig,1).gt.1.E4)then
678	rnat(ig)=0
679	endif
680	enddo
681
682	print*,'WARNING! Surface type currently decided by surface inertia'
683	print*,'This should be improved e.g. in newstart.F'
684	endif!(.not.ok_slab_ocean)
685
686	! initialise surface history variable
687	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
688	qsurf_hist(:,:)=qsurf(:,:)
689
690	! initialise variable for dynamical heating diagnostic
691	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
692	ztprevious(:,:)=pt(:,:)
693
694	! Set temperature just above condensation temperature (for Early Mars)
695	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
696	if(nearco2cond) then
697	write(,)' WARNING! Starting at Tcond+1K'
698	do l=1, nlayer
699	do ig=1,ngrid
700	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
701	-pt(ig,l)) / ptimestep
702	enddo
703	enddo
704	endif
705
706	if(meanOLR)then
707	! to record global radiative balance
708	call system('rm -f rad_bal.out')
709	! to record global mean/max/min temperatures
710	call system('rm -f tem_bal.out')
711	! to record global hydrological balance
712	call system('rm -f h2o_bal.out')
713	endif
714
715	watertest=.false.
716	if(water)then
717	! initialise variables for water cycle
718
719	if(enertest)then
720	watertest = .true.
721	endif
722
723	if(ice_update)then
724	ice_initial(:)=qsurf(:,igcm_h2o_ice)
725	endif
726
727	endif
728	call su_watercycle ! even if we don't have a water cycle, we might
729	! need epsi for the wvp definitions in callcorrk.F
730
731	if (ngrid.ne.1) then ! no need to create a restart file in 1d
732	call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
733	ptimestep,pday+nday,time_phys,area, &
734	albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
735	endif
736
737	endif ! (end of "if firstcall")
738
739	! ---------------------------------------------------
740	! 1.2 Initializations done at every physical timestep:
741	! ---------------------------------------------------
742	! Initialize various variables
743	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
744
745	pdu(1:ngrid,1:nlayer) = 0.0
746	pdv(1:ngrid,1:nlayer) = 0.0
747	if ( .not.nearco2cond ) then
748	pdt(1:ngrid,1:nlayer) = 0.0
749	endif
750	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
751	pdpsrf(1:ngrid) = 0.0
752	zflubid(1:ngrid) = 0.0
753	zdtsurf(1:ngrid) = 0.0
754	dqsurf(1:ngrid,1:nq)= 0.0
755	flux_sens_lat(1:ngrid) = 0.0
756	taux(1:ngrid) = 0.0
757	tauy(1:ngrid) = 0.0
758
759
760
761
762	zday=pday+ptime ! compute time, in sols (and fraction thereof)
763
764	! Compute Stellar Longitude (Ls), and orbital parameters
765	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
766	if (season) then
767	call stellarlong(zday,zls)
768	else
769	call stellarlong(float(day_ini),zls)
770	end if
771
772	call orbite(zls,dist_star,declin,right_ascen)
773	if (tlocked) then
774	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
775	elseif (diurnal) then
776	zlss=-2.pi(zday-.5)
777	else if(diurnal .eqv. .false.) then
778	zlss=9999.
779	endif
780
781
782	! Compute variations of g with latitude (oblate case)
783	!
784	if (oblate .eqv. .false.) then
785	glat(:) = g
786	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
787	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
788
789	call abort
790	else
791
792	gmplanet = MassPlanetgrav1e24
793	do ig=1,ngrid
794	glat(ig)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - lati(ig))))
795	end do
796	endif
797
798	!! write(,) 'lati, glat =',lati(1)180./pi,glat(1), lati(ngrid/2)180./pi, glat(ngrid/2)
799
800
801
802	! Compute geopotential between layers
803	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
804
805	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
806	do l=1,nlayer
807	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
808	enddo
809
810	zzlev(1:ngrid,1)=0.
811	zzlev(1:ngrid,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
812
813	do l=2,nlayer
814	do ig=1,ngrid
815	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
816	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
817	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
818	enddo
819	enddo
820	! Potential temperature calculation may not be the same in physiq and dynamic...
821
822	! Compute potential temperature
823	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
824
825	do l=1,nlayer
826	do ig=1,ngrid
827	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
828	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
829	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
830	massarea(ig,l)=mass(ig,l)*area(ig)
831	enddo
832	enddo
833
834	! Compute vertical velocity (m/s) from vertical mass flux
835	! w = F / (rhoarea) and rho = P/(rT)
836	! but first linearly interpolate mass flux to mid-layers
837	do l=1,nlayer-1
838	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
839	enddo
840	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
841	do l=1,nlayer
842	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
843	(pplay(1:ngrid,l)*area(1:ngrid))
844	enddo
845
846	!-----------------------------------------------------------------------
847	! 2. Compute radiative tendencies
848	!-----------------------------------------------------------------------
849
850	if (callrad) then
851	if( mod(icount-1,iradia).eq.0.or.lastcall) then
852
853	! Compute local stellar zenith angles
854	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
855
856	if (tlocked) then
857	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
858	ztim1=SIN(declin)
859	ztim2=COS(declin)*COS(zlss)
860	ztim3=COS(declin)*SIN(zlss)
861
862	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
863	ztim1,ztim2,ztim3,mu0,fract, flatten)
864
865	elseif (diurnal) then
866	ztim1=SIN(declin)
867	ztim2=COS(declin)COS(2.pi*(zday-.5))
868	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
869
870	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
871	ztim1,ztim2,ztim3,mu0,fract, flatten)
872	else if(diurnal .eqv. .false.) then
873
874	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad,flatten)
875	! WARNING: this function appears not to work in 1D
876
877	endif
878
879	!! Eclipse incoming sunlight (e.g. Saturn ring shadowing)
880
881	if(rings_shadow) then
882	write(,) 'Rings shadow activated'
883	if(diurnal .eqv. .false.) then ! we need to compute the daily average insolation
884	pas = 1./nb_hours
885	ptime_day = 0.
886	fract(:) = 0.
887	ALLOCATE(tmp_fract(ngrid))
888	ALLOCATE(tmp_mu0(ngrid))
889	tmp_fract(:) = 0.
890	eclipse(:) = 0.
891	tmp_mu0(:) = 0.
892
893	do m=1, nb_hours
894	ptime_day = m*pas
895	call stellarlong(pday+ptime_day,tmp_zls)
896	call orbite(tmp_zls,tmp_dist_star,tmp_declin,tmp_right_ascen)
897	ztim1=SIN(tmp_declin)
898	ztim2=COS(tmp_declin)COS(2.pi*(pday+ptime_day-.5))
899	ztim3=-COS(tmp_declin)SIN(2.pi*(pday+ptime_day-.5))
900
901
902	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
903	ztim1,ztim2,ztim3,tmp_mu0,tmp_fract, flatten)
904
905	call rings(ngrid, tmp_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,"Lss","sub solar longitude","deg",0,zlss*180./pi)
2077	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
2078	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
2079	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2080	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2081	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2082	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
2083	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2084	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2085	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2086	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
2087
2088	! Subsurface temperatures
2089	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2090	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
2091
2092	! Oceanic layers
2093	if(ok_slab_ocean) then
2094	call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
2095	call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
2096	call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
2097	call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
2098	call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
2099	call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
2100	call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
2101	call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
2102	call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
2103	call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
2104	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2105	call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
2106	endif! (ok_slab_ocean)
2107
2108	! Total energy balance diagnostics
2109	if(callrad.and.(.not.newtonian))then
2110	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
2111	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2112	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2113	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2114	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
2115	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
2116	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
2117	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
2118	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
2119	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
2120	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
2121	if(ok_slab_ocean) then
2122	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
2123	else
2124	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2125	endif!(ok_slab_ocean)
2126	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2127	endif
2128
2129	if(enertest) then
2130	if (calldifv) then
2131	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
2132	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
2133	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
2134	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
2135	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2136	endif
2137	if (corrk) then
2138	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
2139	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
2140	endif
2141	if(watercond) then
2142	! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
2143	! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
2144	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
2145	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
2146	call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
2147	! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
2148	endif
2149	endif
2150
2151	! Temporary inclusions for heating diagnostics
2152	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2153	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2154	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2155	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2156
2157	! debugging
2158	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2159	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2160
2161	! Output aerosols
2162	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
2163	call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
2164	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
2165	call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
2166	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
2167	call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
2168	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
2169	call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
2170
2171	! Output tracers
2172	if (tracer) then
2173	do iq=1,nq
2174	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2175	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2176	! 'kg m^-2',2,qsurf(1,iq) )
2177	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2178	'kg m^-2',2,qsurf_hist(1,iq) )
2179	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2180	'kg m^-2',2,qcol(1,iq) )
2181
2182	if(watercond.or.CLFvarying)then
2183	call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
2184	call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
2185	call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2186	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2187	call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
2188	endif
2189
2190	if(waterrain)then
2191	call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
2192	call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
2193	endif
2194
2195	if((hydrology).and.(.not.ok_slab_ocean))then
2196	call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
2197	endif
2198
2199	if(ice_update)then
2200	call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
2201	call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
2202	endif
2203
2204	do ig=1,ngrid
2205	if(tau_col(ig).gt.1.e3)then
2206	! print*,'WARNING: tau_col=',tau_col(ig)
2207	! print*,'at ig=',ig,'in PHYSIQ'
2208	endif
2209	end do
2210
2211	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2212
2213	enddo
2214	endif
2215
2216	! output spectrum
2217	if(specOLR.and.corrk)then
2218	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2219	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2220	endif
2221
2222	icount=icount+1
2223
2224	if (lastcall) then
2225
2226	! deallocate gas variables
2227	!$OMP BARRIER
2228	!$OMP MASTER
2229	IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
2230	IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac ) ! both allocated in su_gases.F90
2231	!$OMP END MASTER
2232	!$OMP BARRIER
2233
2234	! deallocate saved arrays
2235	! this is probably not that necessary
2236	! ... but probably a good idea to clean the place before we leave
2237	IF ( ALLOCATED(tsurf)) DEALLOCATE(tsurf)
2238	IF ( ALLOCATED(tsoil)) DEALLOCATE(tsoil)
2239	IF ( ALLOCATED(albedo)) DEALLOCATE(albedo)
2240	IF ( ALLOCATED(albedo0)) DEALLOCATE(albedo0)
2241	IF ( ALLOCATED(rnat)) DEALLOCATE(rnat)
2242	IF ( ALLOCATED(emis)) DEALLOCATE(emis)
2243	IF ( ALLOCATED(dtrad)) DEALLOCATE(dtrad)
2244	IF ( ALLOCATED(fluxrad_sky)) DEALLOCATE(fluxrad_sky)
2245	IF ( ALLOCATED(fluxrad)) DEALLOCATE(fluxrad)
2246	IF ( ALLOCATED(capcal)) DEALLOCATE(capcal)
2247	IF ( ALLOCATED(fluxgrd)) DEALLOCATE(fluxgrd)
2248	IF ( ALLOCATED(qsurf)) DEALLOCATE(qsurf)
2249	IF ( ALLOCATED(q2)) DEALLOCATE(q2)
2250	IF ( ALLOCATED(ztprevious)) DEALLOCATE(ztprevious)
2251	IF ( ALLOCATED(hice)) DEALLOCATE(hice)
2252	IF ( ALLOCATED(cloudfrac)) DEALLOCATE(cloudfrac)
2253	IF ( ALLOCATED(totcloudfrac)) DEALLOCATE(totcloudfrac)
2254	IF ( ALLOCATED(qsurf_hist)) DEALLOCATE(qsurf_hist)
2255	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
2256	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
2257	IF ( ALLOCATED(ice_initial)) DEALLOCATE(ice_initial)
2258	IF ( ALLOCATED(ice_min)) DEALLOCATE(ice_min)
2259
2260	IF ( ALLOCATED(fluxsurf_lw)) DEALLOCATE(fluxsurf_lw)
2261	IF ( ALLOCATED(fluxsurf_sw)) DEALLOCATE(fluxsurf_sw)
2262	IF ( ALLOCATED(fluxtop_lw)) DEALLOCATE(fluxtop_lw)
2263	IF ( ALLOCATED(fluxabs_sw)) DEALLOCATE(fluxabs_sw)
2264	IF ( ALLOCATED(fluxtop_dn)) DEALLOCATE(fluxtop_dn)
2265	IF ( ALLOCATED(fluxdyn)) DEALLOCATE(fluxdyn)
2266	IF ( ALLOCATED(OLR_nu)) DEALLOCATE(OLR_nu)
2267	IF ( ALLOCATED(OSR_nu)) DEALLOCATE(OSR_nu)
2268	IF ( ALLOCATED(sensibFlux)) DEALLOCATE(sensibFlux)
2269	IF ( ALLOCATED(zdtlw)) DEALLOCATE(zdtlw)
2270	IF ( ALLOCATED(zdtsw)) DEALLOCATE(zdtsw)
2271	IF ( ALLOCATED(tau_col)) DEALLOCATE(tau_col)
2272	IF ( ALLOCATED(pctsrf_sic)) DEALLOCATE(pctsrf_sic)
2273	IF ( ALLOCATED(tslab)) DEALLOCATE(tslab)
2274	IF ( ALLOCATED(tsea_ice)) DEALLOCATE(tsea_ice)
2275	IF ( ALLOCATED(sea_ice)) DEALLOCATE(sea_ice)
2276	IF ( ALLOCATED(zmasq)) DEALLOCATE(zmasq)
2277	IF ( ALLOCATED(knindex)) DEALLOCATE(knindex)
2278
2279	!! this is defined in comsaison_h
2280	IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
2281
2282	IF ( ALLOCATED(fract)) DEALLOCATE(fract)
2283
2284
2285	!! this is defined in radcommon_h
2286	IF ( ALLOCATED(eclipse)) DEALLOCATE(eclipse)
2287
2288	!! this is defined in comsoil_h
2289	IF ( ALLOCATED(layer)) DEALLOCATE(layer)
2290	IF ( ALLOCATED(mlayer)) DEALLOCATE(mlayer)
2291	IF ( ALLOCATED(inertiedat)) DEALLOCATE(inertiedat)
2292
2293	!! this is defined in surfdat_h
2294	IF ( ALLOCATED(albedodat)) DEALLOCATE(albedodat)
2295	IF ( ALLOCATED(phisfi)) DEALLOCATE(phisfi)
2296	IF ( ALLOCATED(zmea)) DEALLOCATE(zmea)
2297	IF ( ALLOCATED(zstd)) DEALLOCATE(zstd)
2298	IF ( ALLOCATED(zsig)) DEALLOCATE(zsig)
2299	IF ( ALLOCATED(zgam)) DEALLOCATE(zgam)
2300	IF ( ALLOCATED(zthe)) DEALLOCATE(zthe)
2301	IF ( ALLOCATED(dryness)) DEALLOCATE(dryness)
2302	IF ( ALLOCATED(watercaptag)) DEALLOCATE(watercaptag)
2303
2304	!! this is defined in tracer_h
2305	IF ( ALLOCATED(noms)) DEALLOCATE(noms)
2306	IF ( ALLOCATED(mmol)) DEALLOCATE(mmol)
2307	IF ( ALLOCATED(radius)) DEALLOCATE(radius)
2308	IF ( ALLOCATED(rho_q)) DEALLOCATE(rho_q)
2309	IF ( ALLOCATED(qext)) DEALLOCATE(qext)
2310	IF ( ALLOCATED(alpha_lift)) DEALLOCATE(alpha_lift)
2311	IF ( ALLOCATED(alpha_devil)) DEALLOCATE(alpha_devil)
2312	IF ( ALLOCATED(qextrhor)) DEALLOCATE(qextrhor)
2313	IF ( ALLOCATED(igcm_dustbin)) DEALLOCATE(igcm_dustbin)
2314
2315	!! this is defined in comgeomfi_h
2316	IF ( ALLOCATED(lati)) DEALLOCATE(lati)
2317	IF ( ALLOCATED(long)) DEALLOCATE(long)
2318	IF ( ALLOCATED(area)) DEALLOCATE(area)
2319	IF ( ALLOCATED(sinlat)) DEALLOCATE(sinlat)
2320	IF ( ALLOCATED(coslat)) DEALLOCATE(coslat)
2321	IF ( ALLOCATED(sinlon)) DEALLOCATE(sinlon)
2322	IF ( ALLOCATED(coslon)) DEALLOCATE(coslon)
2323	endif
2324
2325	return
2326	end subroutine physiq

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

Download in other formats:

Original Format