Context Navigation

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

Last change on this file since 298 was 298, checked in by milmd, 10 years ago
Less output messages are written. On 20000 cores it is better. In LMDZ, only master of MPI and OpenMP can write.
Property svn:executable set to ``*
File size: 95.9 KB

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

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

Download in other formats: