Context Navigation

physiq_mod.F @ 1754

Last change on this file since 1754 was 1720, checked in by jaudouard, 8 years ago
Update on CO2 ice clouds scheme
File size: 113.9 KB

Line
1	MODULE physiq_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE physiq(
8	$ ngrid,nlayer,nq
9	$ ,firstcall,lastcall
10	$ ,pday,ptime,ptimestep
11	$ ,pplev,pplay,pphi
12	$ ,pu,pv,pt,pq
13	$ ,flxw
14	$ ,pdu,pdv,pdt,pdq,pdpsrf)
15
16	use watercloud_mod
17	use aeropacity_mod
18	use callradite_mod
19	use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
20	& igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
21	& igcm_ccn_mass, igcm_ccn_number,
22	& igcm_ccnco2_mass, igcm_ccnco2_number,
23	& rho_ice_co2,nuiceco2_sed,nuiceco2_ref,
24	& igcm_dust_mass, igcm_dust_number, igcm_h2o2,
25	& nuice_ref, rho_ice, rho_dust, ref_r0,
26	& igcm_he
27	use comsoil_h, only: inertiedat, ! soil thermal inertia
28	& tsoil, nsoilmx ! number of subsurface layers
29	use geometry_mod, only: longitude, latitude, cell_area
30	use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
31	use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
32	& zthe, z0, albedo_h2o_ice,
33	& frost_albedo_threshold,
34	& tsurf, co2ice, emis,
35	& capcal, fluxgrd, qsurf
36	use comsaison_h, only: dist_sol, declin, mu0, fract
37	use slope_mod, only: theta_sl, psi_sl
38	use conc_mod, only: rnew, cpnew, mmean
39	use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
40	use dimradmars_mod, only: tauscaling, aerosol,
41	& dtrad, fluxrad_sky, fluxrad, albedo,
42	& naerkind
43	use turb_mod, only: q2, wstar, ustar, sensibFlux,
44	& zmax_th, hfmax_th, turb_resolved
45	use planete_h, only: aphelie, periheli, year_day, peri_day,
46	& obliquit
47	USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
48	use param_v4_h, only: nreact,n_avog,
49	& fill_data_thermos, allocate_param_thermos
50	use iono_h, only: allocate_param_iono
51	#ifdef MESOSCALE
52	use comsoil_h, only: mlayer,layer
53	use surfdat_h, only: z0_default
54	use comm_wrf
55	#else
56	use planetwide_mod
57	use phyredem, only: physdem0, physdem1
58	use eofdump_mod, only: eofdump
59	USE vertical_layers_mod, ONLY: ap,bp,aps,bps
60	#endif
61
62	IMPLICIT NONE
63	c=======================================================================
64	c
65	c subject:
66	c --------
67	c
68	c Organisation of the physical parametrisations of the LMD
69	c martian atmospheric general circulation model.
70	c
71	c The GCM can be run without or with tracer transport
72	c depending on the value of Logical "tracer" in file "callphys.def"
73	c Tracers may be water vapor, ice OR chemical species OR dust particles
74	c
75	c SEE comments in initracer.F about numbering of tracer species...
76	c
77	c It includes:
78	c
79	c 1. Initialization:
80	c 1.1 First call initializations
81	c 1.2 Initialization for every call to physiq
82	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
83	c 2. Compute radiative transfer tendencies
84	c (longwave and shortwave) for CO2 and aerosols.
85	c 3. Gravity wave and subgrid scale topography drag :
86	c 4. Vertical diffusion (turbulent mixing):
87	c 5. Convective adjustment
88	c 6. Condensation and sublimation of carbon dioxide.
89	c 7. TRACERS :
90	c 7a. water, water ice, co2 ice (clouds)
91	c 7b. call for photochemistry when tracers are chemical species
92	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
93	c 7d. updates (CO2 pressure variations, surface budget)
94	c 8. Contribution to tendencies due to thermosphere
95	c 9. Surface and sub-surface temperature calculations
96	c 10. Write outputs :
97	c - "startfi", "histfi" (if it's time)
98	c - Saving statistics (if "callstats = .true.")
99	c - Dumping eof (if "calleofdump = .true.")
100	c - Output any needed variables in "diagfi"
101	c 11. Diagnostic: mass conservation of tracers
102	c
103	c author:
104	c -------
105	c Frederic Hourdin 15/10/93
106	c Francois Forget 1994
107	c Christophe Hourdin 02/1997
108	c Subroutine completly rewritten by F.Forget (01/2000)
109	c Introduction of the photochemical module: S. Lebonnois (11/2002)
110	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
111	c Water ice clouds: Franck Montmessin (update 06/2003)
112	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
113	c Nb: See callradite.F for more information.
114	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
115	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
116	c
117	c 10/16 J. Audouard: modifications for CO2 clouds scheme
118
119	c arguments:
120	c ----------
121	c
122	c input:
123	c ------
124	c ecri period (in dynamical timestep) to write output
125	c ngrid Size of the horizontal grid.
126	c All internal loops are performed on that grid.
127	c nlayer Number of vertical layers.
128	c nq Number of advected fields
129	c firstcall True at the first call
130	c lastcall True at the last call
131	c pday Number of days counted from the North. Spring
132	c equinoxe.
133	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
134	c ptimestep timestep (s)
135	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
136	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
137	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
138	c pu(ngrid,nlayer) u component of the wind (ms-1)
139	c pv(ngrid,nlayer) v component of the wind (ms-1)
140	c pt(ngrid,nlayer) Temperature (K)
141	c pq(ngrid,nlayer,nq) Advected fields
142	c pudyn(ngrid,nlayer) \|
143	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
144	c ptdyn(ngrid,nlayer) \| corresponding variables
145	c pqdyn(ngrid,nlayer,nq) \|
146	c flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
147	c
148	c output:
149	c -------
150	c
151	c pdu(ngrid,nlayer) \|
152	c pdv(ngrid,nlayer) \| Temporal derivative of the corresponding
153	c pdt(ngrid,nlayer) \| variables due to physical processes.
154	c pdq(ngrid,nlayer,nq) \|
155	c pdpsrf(ngrid) \|
156
157	c
158	c=======================================================================
159	c
160	c 0. Declarations :
161	c ------------------
162
163	#include "callkeys.h"
164	#include "comg1d.h"
165	#include "nlteparams.h"
166	#include "chimiedata.h"
167	#include "netcdf.inc"
168
169	c Arguments :
170	c -----------
171
172	c inputs:
173	c -------
174	INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
175	INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
176	INTEGER,INTENT(in) :: nq ! number of tracers
177	LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
178	LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
179	REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
180	REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
181	REAL,INTENT(in) :: ptimestep ! physics timestep (s)
182	REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
183	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
184	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
185	REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
186	REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
187	REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
188	REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
189	REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
190	! at lower boundary of layer
191
192	c outputs:
193	c --------
194	c physical tendencies
195	REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
196	REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
197	REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
198	REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
199	REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
200
201
202	c Local saved variables:
203	c ----------------------
204	INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
205	INTEGER,SAVE :: icount ! counter of calls to physiq during the run.
206	#ifdef DUSTSTORM
207	REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
208	#endif
209	c Variables used by the water ice microphysical scheme:
210	REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m)
211	REAL nuice(ngrid,nlayer) ! Estimated effective variance
212	! of the size distribution
213	real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
214	real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3)
215	REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
216	! thermal inertia (J.s-1/2.m-2.K-1)
217	! (used only when tifeedback=.true.)
218	c Variables used by the CO2 clouds microphysical scheme:
219	REAL riceco2(ngrid,nlayer) ! co2 ice geometric mean radius (m)
220	real rsedcloudco2(ngrid,nlayer) !CO2 Cloud sedimentation radius
221	real rhocloudco2(ngrid,nlayer) !co2 Cloud density (kg.m-3)
222	real zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
223	c Variables used by the photochemistry
224	logical :: asis ! true : adaptative semi-implicit symmetric (asis) chemical solver
225	! false : euler backward chemical solver
226	REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
227	REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry)
228	c Variables used by the slope model
229	REAL sl_ls, sl_lct, sl_lat
230	REAL sl_tau, sl_alb, sl_the, sl_psi
231	REAL sl_fl0, sl_flu
232	REAL sl_ra, sl_di0
233	REAL sky
234
235	REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
236
237	c Local variables :
238	c -----------------
239
240	REAL CBRT
241	EXTERNAL CBRT
242
243	! CHARACTER*80 fichier
244	INTEGER l,ig,ierr,igout,iq,tapphys
245
246	REAL fluxsurf_lw(ngrid) !incident LW (IR) surface flux (W.m-2)
247	REAL fluxsurf_sw(ngrid,2) !incident SW (solar) surface flux (W.m-2)
248	REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2)
249	REAL fluxtop_sw(ngrid,2) !Outgoing SW (solar) flux to space (W.m-2)
250	REAL tauref(ngrid) ! Reference column optical depth at odpref
251	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
252	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
253	! AS: TBD: this one should be in a module !
254	REAL dsodust(ngrid,nlayer) ! density-scaled opacity (in infrared)
255	REAL zls ! solar longitude (rad)
256	REAL zday ! date (time since Ls=0, in martian days)
257	REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers
258	REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
259	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
260
261	c Tendancies due to various processes:
262	REAL dqsurf(ngrid,nq)
263	REAL zdtlw(ngrid,nlayer) ! (K/s)
264	REAL zdtsw(ngrid,nlayer) ! (K/s)
265	! REAL cldtlw(ngrid,nlayer) ! (K/s) LW heating rate for clear area
266	! REAL cldtsw(ngrid,nlayer) ! (K/s) SW heating rate for clear area
267	REAL zdtnirco2(ngrid,nlayer) ! (K/s)
268	REAL zdtnlte(ngrid,nlayer) ! (K/s)
269	REAL zdtsurf(ngrid) ! (K/s)
270	REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
271	REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
272	REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid) ! (K/s)
273	REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
274	REAL zdhadj(ngrid,nlayer) ! (K/s)
275	REAL zdtgw(ngrid,nlayer) ! (K/s)
276	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
277	REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
278	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
279
280	REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
281	REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
282	REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
283	REAL zdqadj(ngrid,nlayer,nq)
284	REAL zdqc(ngrid,nlayer,nq)
285	REAL zdqcloud(ngrid,nlayer,nq),zdqcloudco2(ngrid,nlayer,nq)
286	REAL zdqscloud(ngrid,nq)
287	REAL zdqchim(ngrid,nlayer,nq)
288	REAL zdqschim(ngrid,nq)
289
290	REAL zdteuv(ngrid,nlayer) ! (K/s)
291	REAL zdtconduc(ngrid,nlayer) ! (K/s)
292	REAL zdumolvis(ngrid,nlayer)
293	REAL zdvmolvis(ngrid,nlayer)
294	real zdqmoldiff(ngrid,nlayer,nq)
295
296	c Local variable for local intermediate calcul:
297	REAL zflubid(ngrid)
298	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
299	REAL zdum1(ngrid,nlayer)
300	REAL zdum2(ngrid,nlayer)
301	REAL ztim1,ztim2,ztim3, z1,z2
302	REAL ztime_fin
303	REAL zdh(ngrid,nlayer)
304	REAL zh(ngrid,nlayer) ! potential temperature (K)
305	REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
306	INTEGER length
307	PARAMETER (length=100)
308
309	c local variables only used for diagnostic (output in file "diagfi" or "stats")
310	c -----------------------------------------------------------------------------
311	REAL ps(ngrid), zt(ngrid,nlayer)
312	REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
313	REAL zq(ngrid,nlayer,nq)
314	REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
315	character*2 str2
316	! character*5 str5
317	real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
318	real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
319	integer igmin, lmin
320	logical tdiag
321
322	real co2col(ngrid) ! CO2 column
323	! pplev and pplay are dynamical inputs and must not be modified in the physics.
324	! instead, use zplay and zplev :
325	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
326	! REAL zstress(ngrid),cd
327	real tmean, zlocal(nlayer)
328	real rho(ngrid,nlayer) ! density
329	real vmr(ngrid,nlayer) ! volume mixing ratio
330	real rhopart(ngrid,nlayer) ! number density of a given species
331	real colden(ngrid,nq) ! vertical column of tracers
332	real mass(nq) ! global mass of tracers (g)
333	REAL mtot(ngrid) ! Total mass of water vapor (kg/m2)
334	REAL icetot(ngrid) ! Total mass of water ice (kg/m2)
335	REAL mtotco2(ngrid) ! Total mass of co2 vapor (kg/m2)
336	REAL icetotco2(ngrid) ! Total mass of co2 ice (kg/m2)
337	REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2)
338	REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2)
339	REAL rave(ngrid) ! Mean water ice effective radius (m)
340	REAL raveco2(ngrid) ! Mean co2 ice effective radius (m)
341	REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1
342	REAL tauTES(ngrid) ! column optical depth at 825 cm-1
343	REAL Qabsice ! Water ice absorption coefficient
344	REAL taucloudtes(ngrid) ! Cloud opacity at infrared
345	! reference wavelength using
346	! Qabs instead of Qext
347	! (direct comparison with TES)
348
349	REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
350	REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
351	REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
352	REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
353	REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg)
354	REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg)
355	REAL nccnco2(ngrid,nlayer) ! true n ccnco2 (kg/kg)
356	REAL qccnco2(ngrid,nlayer) ! true q ccnco2 (kg/kg)
357
358	c Test 1d/3d scavenging
359	real h2otot(ngrid)
360	REAL satu(ngrid,nlayer) ! satu ratio for output
361	REAL zqsat(ngrid,nlayer) ! saturation
362	REAL satuco2(ngrid,nlayer) ! co2 satu ratio for output
363	REAL zqsatco2(ngrid,nlayer) ! saturation co2
364	REAL,SAVE :: time_phys
365
366	! Added for new NLTE scheme
367
368	real co2vmr_gcm(ngrid,nlayer)
369	real n2vmr_gcm(ngrid,nlayer)
370	real ovmr_gcm(ngrid,nlayer)
371	real covmr_gcm(ngrid,nlayer)
372	integer ierr_nlte
373	real*8 varerr
374
375	c Variables for PBL
376	REAL zz1(ngrid)
377	REAL lmax_th_out(ngrid)
378	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
379	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
380	INTEGER lmax_th(ngrid),dimout,n_out,n
381	CHARACTER(50) zstring
382	REAL dtke_th(ngrid,nlayer+1)
383	REAL zcdv(ngrid), zcdh(ngrid)
384	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
385	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
386	REAL u_out1(ngrid)
387	REAL T_out1(ngrid)
388	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
389	REAL tstar(ngrid) ! friction velocity and friction potential temp
390	REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
391	! REAL zu2(ngrid)
392
393	c sub-grid scale water ice clouds (A. Pottier 2013)
394	logical clearsky
395	! flux for the part without clouds
396	real zdtsw1(ngrid,nlayer)
397	real zdtlw1(ngrid,nlayer)
398	real fluxsurf_lw1(ngrid)
399	real fluxsurf_sw1(ngrid,2)
400	real fluxtop_lw1(ngrid)
401	real fluxtop_sw1(ngrid,2)
402	REAL taucloudtes1(ngrid)
403	! tf: fraction of clouds, ntf: fraction without clouds
404	real tf, ntf
405	real,allocatable,dimension(:),save :: totcloudfrac ! total cloud fraction over the column
406	REAL rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
407
408	c=======================================================================
409
410	c 1. Initialisation:
411	c -----------------
412
413	c 1.1 Initialisation only at first call
414	c ---------------------------------------
415	IF (firstcall) THEN
416	ALLOCATE(totcloudfrac(ngrid)) !A. Pottier: souci avec ca a priori, tentative
417	! de les allouer dans dimradmars_mod.F90
418
419	c variables set to 0
420	c ~~~~~~~~~~~~~~~~~~
421	aerosol(:,:,:)=0
422	dtrad(:,:)=0
423
424	#ifndef MESOSCALE
425	fluxrad(:)=0
426	wstar(:)=0.
427	#endif
428
429	c read startfi
430	c ~~~~~~~~~~~~
431	#ifndef MESOSCALE
432	! GCM. Read netcdf initial physical parameters.
433	CALL phyetat0 ("startfi.nc",0,0,
434	& nsoilmx,ngrid,nlayer,nq,
435	& day_ini,time_phys,
436	& tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling,
437	& totcloudfrac)
438
439	if (pday.ne.day_ini) then
440	write(,) "PHYSIQ: ERROR: bad synchronization between ",
441	& "physics and dynamics"
442	write(,) "dynamics day: ",pday
443	write(,) "physics day: ",day_ini
444	stop
445	endif
446
447	write (,) 'In physiq day_ini =', day_ini
448
449	#else
450	! MESOSCALE. Supposedly everything is already set in modules.
451	! So we just check. And we fill day_ini
452	print*,"check: --- in physiq.F"
453	print*,"check: rad,cpp,g,r,rcp,daysec"
454	print*,rad,cpp,g,r,rcp,daysec
455	PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
456	PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
457	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
458	PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
459	PRINT*,'check: tracernames ', noms
460	PRINT*,'check: emis ',emis(1),emis(ngrid)
461	PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
462	PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
463	PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid)
464	!!!
465	day_ini = pday
466	#endif
467
468	c initialize tracers
469	c ~~~~~~~~~~~~~~~~~~
470	IF (tracer) THEN
471	CALL initracer(ngrid,nq,qsurf)
472	ENDIF ! end tracer
473
474	c Initialize albedo and orbital calculation
475	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
476	CALL surfini(ngrid,co2ice,qsurf,albedo)
477	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
478
479	c initialize soil
480	c ~~~~~~~~~~~~~~~
481	IF (callsoil) THEN
482	c Thermal inertia feedback:
483	IF (tifeedback) THEN
484	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
485	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
486	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
487	ELSE
488	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
489	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
490	ENDIF ! of IF (tifeedback)
491	ELSE
492	PRINT*,
493	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
494	DO ig=1,ngrid
495	capcal(ig)=1.e5
496	fluxgrd(ig)=0.
497	ENDDO
498	ENDIF
499	icount=1
500
501	#ifndef MESOSCALE
502	c Initialize thermospheric parameters
503	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
504
505	if (callthermos) then
506	call fill_data_thermos
507	call allocate_param_thermos(nlayer)
508	call allocate_param_iono(nlayer,nreact)
509	call param_read_e107
510	endif
511	#endif
512	c Initialize R and Cp as constant
513
514	if (.not.callthermos .and. .not.photochem) then
515	do l=1,nlayer
516	do ig=1,ngrid
517	rnew(ig,l)=r
518	cpnew(ig,l)=cpp
519	mmean(ig,l)=mugaz
520	enddo
521	enddo
522	endif
523
524	if(callnlte.and.nltemodel.eq.2) call nlte_setup
525	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
526	if(thermochem) call chemthermos_readini
527
528	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
529	write(,)"physiq: water_param Surface water ice albedo:",
530	. albedo_h2o_ice
531	ENDIF
532
533	#ifndef MESOSCALE
534	if (callslope) call getslopes(ngrid,phisfi)
535
536	if (ngrid.ne.1) then ! no need to create a restart file in 1d
537	call physdem0("restartfi.nc",longitude,latitude,
538	& nsoilmx,ngrid,nlayer,nq,
539	& ptimestep,pday,time_phys,cell_area,
540	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
541	endif
542	#endif
543
544	ENDIF ! (end of "if firstcall")
545
546
547	c ---------------------------------------------------
548	c 1.2 Initializations done at every physical timestep:
549	c ---------------------------------------------------
550	c
551
552	c Initialize various variables
553	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
554	pdv(:,:)=0
555	pdu(:,:)=0
556	pdt(:,:)=0
557	pdq(:,:,:)=0
558	pdpsrf(:)=0
559	zflubid(:)=0
560	zdtsurf(:)=0
561	dqsurf(:,:)=0
562	#ifdef DUSTSTORM
563	pq_tmp(:,:,:)=0
564	#endif
565	igout=ngrid/2+1
566
567
568	zday=pday+ptime ! compute time, in sols (and fraction thereof)
569
570	c Compute Solar Longitude (Ls) :
571	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
572	if (season) then
573	call solarlong(zday,zls)
574	else
575	call solarlong(float(day_ini),zls)
576	end if
577
578	c Initialize pressure levels
579	c ~~~~~~~~~~~~~~~~~~
580	zplev(:,:) = pplev(:,:)
581	zplay(:,:) = pplay(:,:)
582	ps(:) = pplev(:,1)
583
584
585	c Compute geopotential at interlayers
586	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
587	c ponderation des altitudes au niveau des couches en dp/p
588
589	DO l=1,nlayer
590	DO ig=1,ngrid
591	zzlay(ig,l)=pphi(ig,l)/g
592	ENDDO
593	ENDDO
594	DO ig=1,ngrid
595	zzlev(ig,1)=0.
596	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
597	ENDDO
598	DO l=2,nlayer
599	DO ig=1,ngrid
600	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
601	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
602	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
603	ENDDO
604	ENDDO
605
606
607	! Potential temperature calculation not the same in physiq and dynamic
608
609	c Compute potential temperature
610	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
611	DO l=1,nlayer
612	DO ig=1,ngrid
613	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
614	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
615	ENDDO
616	ENDDO
617
618	#ifndef MESOSCALE
619	c-----------------------------------------------------------------------
620	c 1.2.5 Compute mean mass, cp, and R
621	c --------------------------------
622
623	if(photochem.or.callthermos) then
624	call concentrations(ngrid,nlayer,nq,
625	& zplay,pt,pdt,pq,pdq,ptimestep)
626	endif
627	#endif
628
629	! Compute vertical velocity (m/s) from vertical mass flux
630	! w = F / (rhoarea) and rho = P/(rT)
631	! but first linearly interpolate mass flux to mid-layers
632	do l=1,nlayer-1
633	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
634	enddo
635	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
636	do l=1,nlayer
637	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) /
638	& (pplay(1:ngrid,l)*cell_area(1:ngrid))
639	! NB: here we use r and not rnew since this diagnostic comes
640	! from the dynamics
641	enddo
642
643	c-----------------------------------------------------------------------
644	c 2. Compute radiative tendencies :
645	c------------------------------------
646
647
648	IF (callrad) THEN
649	IF( MOD(icount-1,iradia).EQ.0) THEN
650
651	c Local Solar zenith angle
652	c ~~~~~~~~~~~~~~~~~~~~~~~~
653	CALL orbite(zls,dist_sol,declin)
654
655	IF(diurnal) THEN
656	ztim1=SIN(declin)
657	ztim2=COS(declin)COS(2.pi*(zday-.5))
658	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
659
660	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
661	s ztim1,ztim2,ztim3, mu0,fract)
662
663	ELSE
664	CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
665	ENDIF
666
667	c NLTE cooling from CO2 emission
668	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
669	IF(callnlte) then
670	if(nltemodel.eq.0.or.nltemodel.eq.1) then
671	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
672	else if(nltemodel.eq.2) then
673	co2vmr_gcm(1:ngrid,1:nlayer)=
674	& pq(1:ngrid,1:nlayer,igcm_co2)*
675	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
676	n2vmr_gcm(1:ngrid,1:nlayer)=
677	& pq(1:ngrid,1:nlayer,igcm_n2)*
678	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
679	covmr_gcm(1:ngrid,1:nlayer)=
680	& pq(1:ngrid,1:nlayer,igcm_co)*
681	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
682	ovmr_gcm(1:ngrid,1:nlayer)=
683	& pq(1:ngrid,1:nlayer,igcm_o)*
684	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
685
686	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
687	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
688	$ ovmr_gcm, zdtnlte,ierr_nlte,varerr )
689	if(ierr_nlte.gt.0) then
690	write(,)
691	$ 'WARNING: nlte_tcool output with error message',
692	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
693	write(,)'I will continue anyway'
694	endif
695
696	zdtnlte(1:ngrid,1:nlayer)=
697	& zdtnlte(1:ngrid,1:nlayer)/86400.
698	endif
699	else
700	zdtnlte(:,:)=0.
701	endif
702
703	c Find number of layers for LTE radiation calculations
704	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
705	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
706
707	c Note: Dustopacity.F has been transferred to callradite.F
708
709	#ifdef DUSTSTORM
710	!! specific case: save the quantity of dust before adding perturbation
711	if (firstcall) then
712	pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
713	pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
714	endif
715	#endif
716
717	c Call main radiative transfer scheme
718	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
719	c Transfer through CO2 (except NIR CO2 absorption)
720	c and aerosols (dust and water ice)
721
722	c Radiative transfer
723	c ------------------
724	! callradite for the part with clouds
725	clearsky=.false. ! part with clouds for both cases CLFvarying true/false
726	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
727	$ emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
728	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
729	$ tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
730	$ nuice,co2ice,clearsky,totcloudfrac)
731	! case of sub-grid water ice clouds: callradite for the clear case
732	IF (CLFvarying) THEN
733	! ---> PROBLEMS WITH ALLOCATED ARRAYS
734	! (temporary solution in callcorrk: do not deallocate
735	! if
736	! CLFvarying ...) ?? AP ??
737	clearsky=.true. !
738	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
739	& albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
740	& dist_sol,igout,zdtlw1,zdtsw1,fluxsurf_lw1,
741	& fluxsurf_sw1,fluxtop_lw1,fluxtop_sw1,tauref,tau,
742	& aerosol,dsodust,tauscaling,taucloudtes1,rdust,
743	& rice,nuice,co2ice, clearsky, totcloudfrac)
744	clearsky = .false. ! just in case.
745	! Sum the fluxes and heating rates from cloudy/clear
746	! cases
747	DO ig=1,ngrid
748	tf=totcloudfrac(ig)
749	ntf=1.-tf
750	fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig)
751	& + tf*fluxsurf_lw(ig)
752	fluxsurf_sw(ig,1) = ntf*fluxsurf_sw1(ig,1)
753	& + tf*fluxsurf_sw(ig,1)
754	fluxsurf_sw(ig,2) = ntf*fluxsurf_sw1(ig,2)
755	& + tf*fluxsurf_sw(ig,2)
756	fluxtop_lw(ig) = ntf*fluxtop_lw1(ig)
757	& + tf*fluxtop_lw(ig)
758	fluxtop_sw(ig,1) = ntf*fluxtop_sw1(ig,1)
759	& + tf*fluxtop_sw(ig,1)
760	fluxtop_sw(ig,2) = ntf*fluxtop_sw1(ig,2)
761	& + tf*fluxtop_sw(ig,2)
762	taucloudtes(ig) = ntf*taucloudtes1(ig)
763	& + tf*taucloudtes(ig)
764	zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer)
765	& + tf*zdtlw(ig,1:nlayer)
766	zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer)
767	& + tf*zdtsw(ig,1:nlayer)
768	ENDDO
769
770	ENDIF ! (CLFvarying)
771
772	#ifdef DUSTSTORM
773	!! specific case: compute the added quantity of dust for perturbation
774	if (firstcall) then
775	pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
776	& pdq(1:ngrid,1:nlayer,igcm_dust_mass)
777	& - pq_tmp(1:ngrid,1:nlayer,1)
778	& + pq(1:ngrid,1:nlayer,igcm_dust_mass)
779	pdq(1:ngrid,1:nlayer,igcm_dust_number)=
780	& pdq(1:ngrid,1:nlayer,igcm_dust_number)
781	& - pq_tmp(1:ngrid,1:nlayer,2)
782	& + pq(1:ngrid,1:nlayer,igcm_dust_number)
783	endif
784	#endif
785
786	c Outputs for basic check (middle of domain)
787	c ------------------------------------------
788	write(*,'("Ls =",f11.6," check lat =",f10.6,
789	& " lon =",f11.6)')
790	& zls180./pi,latitude(igout)180/pi,
791	& longitude(igout)*180/pi
792	write(*,'(" tauref(",f4.0," Pa) =",f9.6,
793	& " tau(",f4.0," Pa) =",f9.6)')
794	& odpref,tauref(igout),
795	& odpref,tau(igout,1)*odpref/zplev(igout,1)
796	c ---------------------------------------------------------
797	c Call slope parameterization for direct and scattered flux
798	c ---------------------------------------------------------
799	IF(callslope) THEN
800	print *, 'Slope scheme is on and computing...'
801	DO ig=1,ngrid
802	sl_the = theta_sl(ig)
803	IF (sl_the .ne. 0.) THEN
804	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
805	DO l=1,2
806	sl_lct = ptime24. + 180.longitude(ig)/pi/15.
807	sl_ra = pi*(1.0-sl_lct/12.)
808	sl_lat = 180.*latitude(ig)/pi
809	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
810	sl_alb = albedo(ig,l)
811	sl_psi = psi_sl(ig)
812	sl_fl0 = fluxsurf_sw(ig,l)
813	sl_di0 = 0.
814	if (mu0(ig) .gt. 0.) then
815	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
816	sl_di0 = sl_di0*1370./dist_sol/dist_sol
817	sl_di0 = sl_di0/ztim1
818	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
819	endif
820	! you never know (roundup concern...)
821	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
822	!!!!!!!!!!!!!!!!!!!!!!!!!!
823	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
824	& sl_tau, sl_alb, sl_the, sl_psi,
825	& sl_di0, sl_fl0, sl_flu )
826	!!!!!!!!!!!!!!!!!!!!!!!!!!
827	fluxsurf_sw(ig,l) = sl_flu
828	ENDDO
829	!!! compute correction on IR flux as well
830	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
831	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
832	ENDIF
833	ENDDO
834	ENDIF
835
836	c CO2 near infrared absorption
837	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
838	zdtnirco2(:,:)=0
839	if (callnirco2) then
840	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
841	. mu0,fract,declin, zdtnirco2)
842	endif
843
844	c Radiative flux from the sky absorbed by the surface (W.m-2)
845	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
846	DO ig=1,ngrid
847	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
848	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
849	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
850	ENDDO
851
852
853	c Net atmospheric radiative heating rate (K.s-1)
854	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
855	IF(callnlte) THEN
856	CALL blendrad(ngrid, nlayer, zplay,
857	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
858	ELSE
859	DO l=1,nlayer
860	DO ig=1,ngrid
861	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
862	& +zdtnirco2(ig,l)
863	ENDDO
864	ENDDO
865	ENDIF
866
867	ENDIF ! of if(mod(icount-1,iradia).eq.0)
868
869	c Transformation of the radiative tendencies:
870	c -------------------------------------------
871
872	c Net radiative surface flux (W.m-2)
873	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
874	c
875	DO ig=1,ngrid
876	zplanck(ig)=tsurf(ig)*tsurf(ig)
877	zplanck(ig)=emis(ig)*
878	$ stephanzplanck(ig)zplanck(ig)
879	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
880	IF(callslope) THEN
881	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
882	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
883	ENDIF
884	ENDDO
885
886	DO l=1,nlayer
887	DO ig=1,ngrid
888	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
889	ENDDO
890	ENDDO
891
892	ENDIF ! of IF (callrad)
893
894	c-----------------------------------------------------------------------
895	c 3. Gravity wave and subgrid scale topography drag :
896	c -------------------------------------------------
897
898
899	IF(calllott)THEN
900
901	CALL calldrag_noro(ngrid,nlayer,ptimestep,
902	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
903
904	DO l=1,nlayer
905	DO ig=1,ngrid
906	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
907	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
908	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
909	ENDDO
910	ENDDO
911	ENDIF
912
913	c-----------------------------------------------------------------------
914	c 4. Vertical diffusion (turbulent mixing):
915	c -----------------------------------------
916
917	IF (calldifv) THEN
918
919	DO ig=1,ngrid
920	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
921	ENDDO
922
923	zdum1(:,:)=0
924	zdum2(:,:)=0
925	do l=1,nlayer
926	do ig=1,ngrid
927	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
928	enddo
929	enddo
930
931	c ----------------------
932	c Treatment of a special case : using new surface layer (Richardson based)
933	c without using the thermals in gcm and mesoscale can yield problems in
934	c weakly unstable situations when winds are near to 0. For those cases, we add
935	c a unit subgrid gustiness. Remember that thermals should be used we using the
936	c Richardson based surface layer model.
937	IF ( .not.calltherm
938	. .and. callrichsl
939	. .and. .not.turb_resolved) THEN
940	DO ig=1, ngrid
941	IF (zh(ig,1) .lt. tsurf(ig)) THEN
942	wstar(ig)=1.
943	hfmax_th(ig)=0.2
944	ELSE
945	wstar(ig)=0.
946	hfmax_th(ig)=0.
947	ENDIF
948	ENDDO
949	ENDIF
950	c ----------------------
951
952	IF (tke_heat_flux .ne. 0.) THEN
953	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
954	& (-alog(zplay(:,1)/zplev(:,1)))
955	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
956	ENDIF
957
958	c Calling vdif (Martian version WITH CO2 condensation)
959	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
960	$ ptimestep,capcal,lwrite,
961	$ zplay,zplev,zzlay,zzlev,z0,
962	$ pu,pv,zh,pq,tsurf,emis,qsurf,
963	$ zdum1,zdum2,zdh,pdq,zflubid,
964	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
965	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux)
966
967
968	DO ig=1,ngrid
969	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
970	ENDDO
971
972	IF (.not.turb_resolved) THEN
973	DO l=1,nlayer
974	DO ig=1,ngrid
975	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
976	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
977	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
978
979	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
980	ENDDO
981	ENDDO
982
983	if (tracer) then
984	DO iq=1, nq
985	DO l=1,nlayer
986	DO ig=1,ngrid
987	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
988	ENDDO
989	ENDDO
990	ENDDO
991	DO iq=1, nq
992	DO ig=1,ngrid
993	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
994	ENDDO
995	ENDDO
996	end if ! of if (tracer)
997	ELSE
998	write (,) '******************************************'
999	write (,) '** LES mode: the difv part is only used to'
1000	write (,) '** - provide HFX and UST to the dynamics'
1001	write (,) '** - update TSURF'
1002	write (,) '******************************************'
1003	!! Specific treatment for lifting in turbulent-resolving mode (AC)
1004	IF (lifting .and. doubleq) THEN
1005	!! lifted dust is injected in the first layer.
1006	!! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
1007	!! => lifted dust is not incremented before the sedimentation step.
1008	zdqdif(1:ngrid,1,1:nq)=0.
1009	zdqdif(1:ngrid,1,igcm_dust_number) =
1010	. -zdqsdif(1:ngrid,igcm_dust_number)
1011	zdqdif(1:ngrid,1,igcm_dust_mass) =
1012	. -zdqsdif(1:ngrid,igcm_dust_mass)
1013	zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
1014	DO iq=1, nq
1015	IF ((iq .ne. igcm_dust_mass)
1016	& .and. (iq .ne. igcm_dust_number)) THEN
1017	zdqsdif(:,iq)=0.
1018	ENDIF
1019	ENDDO
1020	ELSE
1021	zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
1022	zdqsdif(1:ngrid,1:nq) = 0.
1023	ENDIF
1024	ENDIF
1025	ELSE
1026	DO ig=1,ngrid
1027	zdtsurf(ig)=zdtsurf(ig)+
1028	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
1029	ENDDO
1030	IF (turb_resolved) THEN
1031	write(,) 'Turbulent-resolving mode !'
1032	write(,) 'Please set calldifv to T in callphys.def'
1033	STOP
1034	ENDIF
1035	ENDIF ! of IF (calldifv)
1036
1037	c-----------------------------------------------------------------------
1038	c 5. Thermals :
1039	c -----------------------------
1040
1041	if(calltherm .and. .not.turb_resolved) then
1042
1043	call calltherm_interface(ngrid,nlayer,nq,
1044	$ tracer,igcm_co2,
1045	$ zzlev,zzlay,
1046	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
1047	$ zplay,zplev,pphi,zpopsk,
1048	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
1049	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
1050
1051	DO l=1,nlayer
1052	DO ig=1,ngrid
1053	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
1054	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
1055	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
1056	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
1057	ENDDO
1058	ENDDO
1059
1060	DO ig=1,ngrid
1061	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
1062	ENDDO
1063
1064	if (tracer) then
1065	DO iq=1,nq
1066	DO l=1,nlayer
1067	DO ig=1,ngrid
1068	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
1069	ENDDO
1070	ENDDO
1071	ENDDO
1072	endif
1073
1074	lmax_th_out(:)=real(lmax_th(:))
1075
1076	else !of if calltherm
1077	lmax_th(:)=0
1078	wstar(:)=0.
1079	hfmax_th(:)=0.
1080	lmax_th_out(:)=0.
1081	end if
1082
1083	c-----------------------------------------------------------------------
1084	c 5. Dry convective adjustment:
1085	c -----------------------------
1086
1087	IF(calladj) THEN
1088
1089	DO l=1,nlayer
1090	DO ig=1,ngrid
1091	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1092	ENDDO
1093	ENDDO
1094	zduadj(:,:)=0
1095	zdvadj(:,:)=0
1096	zdhadj(:,:)=0
1097
1098	CALL convadj(ngrid,nlayer,nq,ptimestep,
1099	$ zplay,zplev,zpopsk,lmax_th,
1100	$ pu,pv,zh,pq,
1101	$ pdu,pdv,zdh,pdq,
1102	$ zduadj,zdvadj,zdhadj,
1103	$ zdqadj)
1104
1105
1106	DO l=1,nlayer
1107	DO ig=1,ngrid
1108	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
1109	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
1110	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
1111
1112	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1113	ENDDO
1114	ENDDO
1115
1116	if(tracer) then
1117	DO iq=1, nq
1118	DO l=1,nlayer
1119	DO ig=1,ngrid
1120	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
1121	ENDDO
1122	ENDDO
1123	ENDDO
1124	end if
1125	ENDIF ! of IF(calladj)
1126
1127
1128
1129	c-----------------------------------------------------------------------
1130	c 6. Specific parameterizations for tracers
1131	c: -----------------------------------------
1132
1133	if (tracer) then
1134
1135	c 6a. Water and ice
1136	c ---------------
1137
1138	c ---------------------------------------
1139	c Water ice condensation in the atmosphere
1140	c ----------------------------------------
1141	IF (water) THEN
1142
1143	call watercloud(ngrid,nlayer,ptimestep,
1144	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
1145	& pq,pdq,zdqcloud,zdtcloud,
1146	& nq,tau,tauscaling,rdust,rice,nuice,
1147	& rsedcloud,rhocloud,totcloudfrac)
1148
1149	c Temperature variation due to latent heat release
1150	if (activice) then
1151	pdt(1:ngrid,1:nlayer) =
1152	& pdt(1:ngrid,1:nlayer) +
1153	& zdtcloud(1:ngrid,1:nlayer)
1154	endif
1155
1156	! increment water vapour and ice atmospheric tracers tendencies
1157	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1158	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1159	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1160	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1161	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1162	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1163
1164	! increment dust and ccn masses and numbers
1165	! We need to check that we have Nccn & Ndust > 0
1166	! This is due to single precision rounding problems
1167	if (microphys) then
1168	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1169	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1170	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1171	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1172	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1173	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1174	where (pq(:,:,igcm_ccn_mass) +
1175	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1176	pdq(:,:,igcm_ccn_mass) =
1177	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1178	pdq(:,:,igcm_ccn_number) =
1179	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1180	end where
1181	where (pq(:,:,igcm_ccn_number) +
1182	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1183	pdq(:,:,igcm_ccn_mass) =
1184	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1185	pdq(:,:,igcm_ccn_number) =
1186	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1187	end where
1188	endif
1189
1190	if (scavenging) then
1191	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1192	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1193	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1194	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1195	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1196	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1197	where (pq(:,:,igcm_dust_mass) +
1198	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1199	pdq(:,:,igcm_dust_mass) =
1200	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1201	pdq(:,:,igcm_dust_number) =
1202	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1203	end where
1204	where (pq(:,:,igcm_dust_number) +
1205	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1206	pdq(:,:,igcm_dust_mass) =
1207	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1208	pdq(:,:,igcm_dust_number) =
1209	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1210	end where
1211	endif ! of if scavenging
1212
1213
1214	END IF ! of IF (water)
1215
1216	c 6a bis. CO2 clouds (CL & JA)
1217
1218	c ---------------------------------------
1219	c CO2 ice cloud condensation in the atmosphere
1220	c ----------------------------------------
1221
1222
1223	IF (co2clouds) THEN
1224
1225	call co2cloud(ngrid,nlayer,ptimestep,
1226	& zplev,zplay,pdpsrf,zzlay,pt,pdt,
1227	& pq,pdq,zdqcloudco2,zdtcloudco2,
1228	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
1229	& rsedcloudco2,rhocloudco2,
1230	& rsedcloud,rhocloud,zzlev,zdqssed_co2)
1231
1232
1233	c Temperature variation due to latent heat release
1234	if (activice) then !Maybe create activice_co2 ?
1235	pdt(1:ngrid,1:nlayer) =
1236	& pdt(1:ngrid,1:nlayer) +
1237	& zdtcloudco2(1:ngrid,1:nlayer)
1238	endif
1239
1240
1241	! increment dust and ccn masses and numbers
1242	! We need to check that we have Nccn & Ndust > 0
1243	! This is due to single precision rounding problems
1244
1245	pdq(1:ngrid,1:nlayer,igcm_co2) =
1246	& pdq(1:ngrid,1:nlayer,igcm_co2) +
1247	& zdqcloudco2(1:ngrid,1:nlayer,igcm_co2)
1248	pdq(1:ngrid,1:nlayer,igcm_co2_ice) =
1249	& pdq(1:ngrid,1:nlayer,igcm_co2_ice) +
1250	& zdqcloudco2(1:ngrid,1:nlayer,igcm_co2_ice)
1251	pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) =
1252	& pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) +
1253	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_mass)
1254	pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) =
1255	& pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) +
1256	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number)
1257	!Update water ice clouds values as well
1258	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1259	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1260	& zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
1261
1262	! increment dust and ccn masses and numbers
1263	! We need to check that we have Nccn & Ndust > 0
1264	! This is due to single precision rounding problems
1265	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1266	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1267	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
1268	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1269	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1270	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
1271	where (pq(:,:,igcm_ccn_mass) +
1272	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1273	pdq(:,:,igcm_ccn_mass) =
1274	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1275	pdq(:,:,igcm_ccn_number) =
1276	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1277	end where
1278	where (pq(:,:,igcm_ccn_number) +
1279	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1280	pdq(:,:,igcm_ccn_mass) =
1281	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1282	pdq(:,:,igcm_ccn_number) =
1283	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1284	end where
1285	c Negative values?
1286	where (pq(:,:,igcm_ccnco2_mass) +
1287	& ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
1288	pdq(:,:,igcm_ccnco2_mass) =
1289	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1290	pdq(:,:,igcm_ccnco2_number) =
1291	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1292	end where
1293	where (pq(:,:,igcm_ccnco2_number) +
1294	& ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
1295	pdq(:,:,igcm_ccnco2_mass) =
1296	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1297	pdq(:,:,igcm_ccnco2_number) =
1298	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1299	end where
1300
1301	! increment dust tracers tendancies
1302	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1303	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1304	& zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_mass)
1305	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1306	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1307	& zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_number)
1308	c Negative values?
1309	where (pq(:,:,igcm_dust_mass) +
1310	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1311	pdq(:,:,igcm_dust_mass) =
1312	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1313	pdq(:,:,igcm_dust_number) =
1314	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1315	end where
1316	where (pq(:,:,igcm_dust_number) +
1317	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1318	pdq(:,:,igcm_dust_mass) =
1319	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1320	pdq(:,:,igcm_dust_number) =
1321	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1322	end where
1323
1324	END IF ! of IF (co2clouds)
1325
1326
1327	c 6b. Aerosol particles
1328	c -------------------
1329
1330	c ----------
1331	c Dust devil :
1332	c ----------
1333	IF(callddevil) then
1334	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
1335	& zdqdev,zdqsdev)
1336
1337	if (dustbin.ge.1) then
1338	do iq=1,nq
1339	DO l=1,nlayer
1340	DO ig=1,ngrid
1341	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1342	ENDDO
1343	ENDDO
1344	enddo
1345	do iq=1,nq
1346	DO ig=1,ngrid
1347	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1348	ENDDO
1349	enddo
1350	endif ! of if (dustbin.ge.1)
1351
1352	END IF ! of IF (callddevil)
1353
1354	c -------------
1355	c Sedimentation : acts also on water ice
1356	c -------------
1357	IF (sedimentation) THEN
1358	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1359	zdqssed(1:ngrid,1:nq)=0
1360
1361	cSedimentation for co2 clouds tracers are inside co2cloud microtimestep
1362	cZdqssed isn't
1363	call callsedim(ngrid,nlayer, ptimestep,
1364	& zplev,zzlev, zzlay, pt, pdt, rdust, rice,
1365	& rsedcloud,rhocloud,
1366	& pq, pdq, zdqsed, zdqssed,nq,
1367	& tau,tauscaling)
1368	!Flux at the surface of co2 ice computed in co2cloud microtimestep
1369	IF (co2clouds) THEN
1370	zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
1371	ENDIF
1372
1373	DO iq=1, nq
1374	DO l=1,nlayer
1375	DO ig=1,ngrid
1376	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1377	ENDDO
1378	ENDDO
1379	ENDDO
1380	DO iq=1, nq
1381	DO ig=1,ngrid
1382	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1383	ENDDO
1384	ENDDO
1385	END IF ! of IF (sedimentation)
1386
1387	c Add lifted dust to tendancies after sedimentation in the LES (AC)
1388	IF (turb_resolved) THEN
1389	DO iq=1, nq
1390	DO l=1,nlayer
1391	DO ig=1,ngrid
1392	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1393	ENDDO
1394	ENDDO
1395	ENDDO
1396	DO iq=1, nq
1397	DO ig=1,ngrid
1398	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
1399	ENDDO
1400	ENDDO
1401	ENDIF
1402
1403	c
1404	c 6c. Chemical species
1405	c ------------------
1406
1407	#ifndef MESOSCALE
1408	c --------------
1409	c photochemistry :
1410	c --------------
1411	IF (photochem .or. thermochem) then
1412
1413	! dust and ice surface area
1414	call surfacearea(ngrid, nlayer, naerkind,
1415	$ ptimestep, zplay, zzlay,
1416	$ pt, pq, pdq, nq,
1417	$ rdust, rice, tau, tauscaling,
1418	$ surfdust, surfice)
1419	! call photochemistry
1420	! asis = .false. euler-backward scheme
1421	! asis = .true. asis scheme
1422
1423	asis = .true.
1424
1425	if (asis) then
1426	call calchim_asis(ngrid,nlayer,nq,
1427	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
1428	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
1429	$ zdqcloud,zdqscloud,tauref,co2ice,
1430	$ pu,pdu,pv,pdv,surfdust,surfice)
1431	else
1432	call calchim(ngrid,nlayer,nq,
1433	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
1434	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
1435	$ zdqcloud,zdqscloud,tauref,co2ice,
1436	$ pu,pdu,pv,pdv,surfdust,surfice)
1437	end if
1438
1439	! increment values of tracers:
1440	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1441	! tracers is zero anyways
1442	DO l=1,nlayer
1443	DO ig=1,ngrid
1444	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1445	ENDDO
1446	ENDDO
1447	ENDDO ! of DO iq=1,nq
1448
1449	! add condensation tendency for H2O2
1450	if (igcm_h2o2.ne.0) then
1451	DO l=1,nlayer
1452	DO ig=1,ngrid
1453	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1454	& +zdqcloud(ig,l,igcm_h2o2)
1455	ENDDO
1456	ENDDO
1457	endif
1458
1459	! increment surface values of tracers:
1460	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1461	! tracers is zero anyways
1462	DO ig=1,ngrid
1463	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1464	ENDDO
1465	ENDDO ! of DO iq=1,nq
1466
1467	! add condensation tendency for H2O2
1468	if (igcm_h2o2.ne.0) then
1469	DO ig=1,ngrid
1470	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1471	& +zdqscloud(ig,igcm_h2o2)
1472	ENDDO
1473	endif
1474
1475	END IF ! of IF (photochem.or.thermochem)
1476	#endif
1477
1478	c 6d. Updates
1479	c ---------
1480
1481	DO iq=1, nq
1482	DO ig=1,ngrid
1483
1484	c ---------------------------------
1485	c Updating tracer budget on surface
1486	c ---------------------------------
1487	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1488
1489	ENDDO ! (ig)
1490	ENDDO ! (iq)
1491
1492	endif ! of if (tracer)
1493
1494	#ifndef MESOSCALE
1495	c-----------------------------------------------------------------------
1496	c 7. THERMOSPHERE CALCULATION
1497	c-----------------------------------------------------------------------
1498
1499	if (callthermos) then
1500	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
1501	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1502	& pt,pq,pu,pv,pdt,pdq,
1503	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1504
1505	DO l=1,nlayer
1506	DO ig=1,ngrid
1507	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1508	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1509	& +zdteuv(ig,l)
1510	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1511	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1512	DO iq=1, nq
1513	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1514	ENDDO
1515	ENDDO
1516	ENDDO
1517
1518	endif ! of if (callthermos)
1519	#endif
1520
1521	c-----------------------------------------------------------------------
1522	c 8. Carbon dioxide condensation-sublimation:
1523	c (should be the last atmospherical physical process to be computed)
1524	c -------------------------------------------
1525
1526	IF (tituscap) THEN
1527	!!! get the actual co2 seasonal cap from Titus observations
1528	CALL geticecover( ngrid, 180.*zls/pi,
1529	. 180.longitude/pi, 180.latitude/pi, co2ice )
1530	co2ice = co2ice * 10000.
1531	ENDIF
1532
1533
1534	pdpsrf(:) = 0
1535
1536	IF (callcond) THEN
1537	CALL newcondens(ngrid,nlayer,nq,ptimestep,
1538	$ capcal,zplay,zplev,tsurf,pt,
1539	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
1540	$ co2ice,albedo,emis,
1541	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
1542	$ fluxsurf_sw,zls)
1543
1544	DO l=1,nlayer
1545	DO ig=1,ngrid
1546	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1547	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1548	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1549	ENDDO
1550	ENDDO
1551	DO ig=1,ngrid
1552	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1553	ENDDO
1554
1555	IF (tracer) THEN
1556	DO iq=1, nq
1557	DO l=1,nlayer
1558	DO ig=1,ngrid
1559	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1560	ENDDO
1561	ENDDO
1562	ENDDO
1563	ENDIF ! of IF (tracer)
1564
1565	#ifndef MESOSCALE
1566	! update surface pressure
1567	DO ig=1,ngrid
1568	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
1569	ENDDO
1570
1571	! update pressure levels
1572	DO l=1,nlayer
1573	DO ig=1,ngrid
1574	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
1575	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
1576	ENDDO
1577	ENDDO
1578	zplev(:,nlayer+1) = 0.
1579
1580	! update layers altitude
1581	DO l=2,nlayer
1582	DO ig=1,ngrid
1583	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
1584	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
1585	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
1586	ENDDO
1587	ENDDO
1588	#endif
1589
1590	ENDIF ! of IF (callcond)
1591
1592
1593	c-----------------------------------------------------------------------
1594	c 9. Surface and sub-surface soil temperature
1595	c-----------------------------------------------------------------------
1596	c
1597	c
1598	c 9.1 Increment Surface temperature:
1599	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1600
1601	DO ig=1,ngrid
1602	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1603	ENDDO
1604
1605	c Prescribe a cold trap at south pole (except at high obliquity !!)
1606	c Temperature at the surface is set there to be the temperature
1607	c corresponding to equilibrium temperature between phases of CO2
1608
1609
1610	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
1611	!#ifndef MESOSCALE
1612	! if (caps.and.(obliquit.lt.27.)) then => now done in newcondens
1613	! NB: Updated surface pressure, at grid point 'ngrid', is
1614	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1615	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1616	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1617	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
1618	! endif
1619	!#endif
1620	c -------------------------------------------------------------
1621	c Change of surface albedo in case of ground frost
1622	c everywhere except on the north permanent cap and in regions
1623	c covered by dry ice.
1624	c ALWAYS PLACE these lines after newcondens !!!
1625	c -------------------------------------------------------------
1626	do ig=1,ngrid
1627	if ((co2ice(ig).eq.0).and.
1628	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
1629	albedo(ig,1) = albedo_h2o_ice
1630	albedo(ig,2) = albedo_h2o_ice
1631	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
1632	c write(,) "physiq.F frost :"
1633	c & ,latitude(ig)180./pi, longitude(ig)180./pi
1634	endif
1635	enddo ! of do ig=1,ngrid
1636	ENDIF ! of IF (tracer.AND.water.AND.(ngrid.NE.1))
1637
1638	c
1639	c 9.2 Compute soil temperatures and subsurface heat flux:
1640	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1641	IF (callsoil) THEN
1642	c Thermal inertia feedback
1643	IF (tifeedback) THEN
1644	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
1645	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
1646	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
1647	ELSE
1648	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1649	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
1650	ENDIF
1651	ENDIF
1652
1653
1654	c-----------------------------------------------------------------------
1655	c 10. Write output files
1656	c ----------------------
1657
1658	c -------------------------------
1659	c Dynamical fields incrementation
1660	c -------------------------------
1661	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1662	! temperature, zonal and meridional wind
1663	DO l=1,nlayer
1664	DO ig=1,ngrid
1665	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1666	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1667	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1668	ENDDO
1669	ENDDO
1670
1671	! tracers
1672	DO iq=1, nq
1673	DO l=1,nlayer
1674	DO ig=1,ngrid
1675	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1676	ENDDO
1677	ENDDO
1678	ENDDO
1679
1680	! Density
1681	DO l=1,nlayer
1682	DO ig=1,ngrid
1683	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1684	ENDDO
1685	ENDDO
1686
1687	! Potential Temperature
1688
1689	DO ig=1,ngrid
1690	DO l=1,nlayer
1691	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
1692	ENDDO
1693	ENDDO
1694
1695
1696	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1697	c DO ig=1,ngrid
1698	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1699	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1700	c ENDDO
1701
1702	c Sum of fluxes in solar spectral bands (for output only)
1703	DO ig=1,ngrid
1704	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1705	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1706	ENDDO
1707	c ***** TEST ****************************************************
1708	ztim1 = 999
1709	DO l=1,nlayer
1710	DO ig=1,ngrid
1711	if (pt(ig,l).lt.ztim1) then
1712	ztim1 = pt(ig,l)
1713	igmin = ig
1714	lmin = l
1715	end if
1716	ENDDO
1717	ENDDO
1718	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1719	write(,) 'PHYSIQ: stability WARNING :'
1720	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1721	& 'ig l =', igmin, lmin
1722	end if
1723	c *******************************************************************
1724
1725	c ---------------------
1726	c Outputs to the screen
1727	c ---------------------
1728
1729	IF (lwrite) THEN
1730	PRINT*,'Global diagnostics for the physics'
1731	PRINT,'Variables and their increments x and dx/dt dt'
1732	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1733	WRITE(*,'(2f10.5,2f15.5)')
1734	s tsurf(igout),zdtsurf(igout)*ptimestep,
1735	s zplev(igout,1),pdpsrf(igout)*ptimestep
1736	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1737	WRITE(*,'(i4,6f10.5)') (l,
1738	s pu(igout,l),pdu(igout,l)*ptimestep,
1739	s pv(igout,l),pdv(igout,l)*ptimestep,
1740	s pt(igout,l),pdt(igout,l)*ptimestep,
1741	s l=1,nlayer)
1742	ENDIF ! of IF (lwrite)
1743
1744	c ----------------------------------------------------------
1745	c ----------------------------------------------------------
1746	c INTERPOLATIONS IN THE SURFACE-LAYER
1747	c ----------------------------------------------------------
1748	c ----------------------------------------------------------
1749
1750	n_out=0 ! number of elements in the z_out array.
1751	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
1752	! to 5
1753	IF (n_out .ne. 0) THEN
1754
1755	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
1756	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
1757	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
1758
1759	z_out(:)=[3.,2.,1.,0.5,0.1]
1760	u_out(:,:)=0.
1761	T_out(:,:)=0.
1762
1763	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
1764	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
1765	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
1766	! pourquoi ustar recalcule ici? fait dans vdifc.
1767
1768	#ifndef MESOSCALE
1769	IF (ngrid .eq. 1) THEN
1770	dimout=0
1771	ELSE
1772	dimout=2
1773	ENDIF
1774	DO n=1,n_out
1775	write(zstring, '(F8.6)') z_out(n)
1776	call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
1777	& 'potential temperature at z_out','K',dimout,T_out(:,n))
1778	call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
1779	& 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
1780	ENDDO
1781	call WRITEDIAGFI(ngrid,'u_star',
1782	& 'friction velocity','m/s',dimout,ustar)
1783	call WRITEDIAGFI(ngrid,'teta_star',
1784	& 'friction potential temperature','K',dimout,tstar)
1785	! call WRITEDIAGFI(ngrid,'L',
1786	! & 'Monin Obukhov length','m',dimout,L_mo)
1787	call WRITEDIAGFI(ngrid,'vvv',
1788	& 'Vertical velocity variance at zout','m',dimout,vvv)
1789	call WRITEDIAGFI(ngrid,'vhf',
1790	& 'Vertical heat flux at zout','m',dimout,vhf)
1791	#else
1792	T_out1(:)=T_out(:,1)
1793	u_out1(:)=u_out(:,1)
1794	#endif
1795
1796	ENDIF
1797
1798	c ----------------------------------------------------------
1799	c ----------------------------------------------------------
1800	c END OF SURFACE LAYER INTERPOLATIONS
1801	c ----------------------------------------------------------
1802	c ----------------------------------------------------------
1803
1804	IF (ngrid.NE.1) THEN
1805
1806	#ifndef MESOSCALE
1807	c -------------------------------------------------------------------
1808	c Writing NetCDF file "RESTARTFI" at the end of the run
1809	c -------------------------------------------------------------------
1810	c Note: 'restartfi' is stored just before dynamics are stored
1811	c in 'restart'. Between now and the writting of 'restart',
1812	c there will have been the itau=itau+1 instruction and
1813	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1814	c thus we store for time=time+dtvr
1815
1816	IF( ((ecritstart.GT.0) .and.
1817	. (MOD(icount*iphysiq,ecritstart).EQ.0))
1818	. .or. lastcall ) THEN
1819
1820	ztime_fin = pday + ptime + ptimestep/(float(iphysiq)*daysec)
1821	. - day_ini - time_phys
1822	print*, pday,ptime,day_ini, time_phys
1823	write(*,'(A,I7,A,F12.5)')
1824	. 'PHYSIQ: Ecriture du fichier restartfi ; icount=',
1825	. icount,' date=',ztime_fin
1826
1827
1828	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
1829	. ptimestep,ztime_fin,
1830	. tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
1831	. totcloudfrac)
1832
1833	ENDIF
1834	#endif
1835
1836	c -------------------------------------------------------------------
1837	c Calculation of diagnostic variables written in both stats and
1838	c diagfi files
1839	c -------------------------------------------------------------------
1840
1841	if (tracer) then
1842
1843	if(doubleq) then
1844	do ig=1,ngrid
1845	dqdustsurf(ig) =
1846	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
1847	dndustsurf(ig) =
1848	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
1849	ndust(ig,:) =
1850	& zq(ig,:,igcm_dust_number)*tauscaling(ig)
1851	qdust(ig,:) =
1852	& zq(ig,:,igcm_dust_mass)*tauscaling(ig)
1853	enddo
1854	if (scavenging) then
1855	do ig=1,ngrid
1856	dqdustsurf(ig) = dqdustsurf(ig) +
1857	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
1858	dndustsurf(ig) = dndustsurf(ig) +
1859	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
1860	nccn(ig,:) =
1861	& zq(ig,:,igcm_ccn_number)*tauscaling(ig)
1862	qccn(ig,:) =
1863	& zq(ig,:,igcm_ccn_mass)*tauscaling(ig)
1864	enddo
1865	endif
1866	endif
1867	if (co2clouds) then
1868	do ig=1,ngrid
1869	nccnco2(ig,:) =
1870	& zq(ig,:,igcm_ccnco2_number)*tauscaling(ig)
1871	qccnco2(ig,:) =
1872	& zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
1873	enddo
1874
1875
1876	call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
1877	& 'kg/kg',3,qccnco2)
1878	call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
1879	& 'part/kg',3,nccnco2)
1880	call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
1881	& 3,pq(:,:,igcm_co2_ice))
1882	call WRITEDIAGFI(ngrid,'co2','co2','kg/kg',
1883	& 3,pq(:,:,igcm_co2))
1884	call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
1885	& 3,pq(:,:,igcm_h2o_ice))
1886	endif
1887
1888
1889	if (water) then
1890	mtot(:)=0
1891	icetot(:)=0
1892	rave(:)=0
1893	tauTES(:)=0
1894	do ig=1,ngrid
1895	do l=1,nlayer
1896	mtot(ig) = mtot(ig) +
1897	& zq(ig,l,igcm_h2o_vap) *
1898	& (zplev(ig,l) - zplev(ig,l+1)) / g
1899	icetot(ig) = icetot(ig) +
1900	& zq(ig,l,igcm_h2o_ice) *
1901	& (zplev(ig,l) - zplev(ig,l+1)) / g
1902	c Computing abs optical depth at 825 cm-1 in each
1903	c layer to simulate NEW TES retrieval
1904	Qabsice = min(
1905	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1906	& )
1907	opTES(ig,l)= 0.75 * Qabsice *
1908	& zq(ig,l,igcm_h2o_ice) *
1909	& (zplev(ig,l) - zplev(ig,l+1)) / g
1910	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1911	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1912	enddo
1913	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
1914	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1915	enddo
1916	call watersat(ngrid*nlayer,zt,zplay,zqsat)
1917	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
1918
1919	if (scavenging) then
1920	Nccntot(:)= 0
1921	Mccntot(:)= 0
1922	rave(:)=0
1923	do ig=1,ngrid
1924	do l=1,nlayer
1925	Nccntot(ig) = Nccntot(ig) +
1926	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
1927	& *(zplev(ig,l) - zplev(ig,l+1)) / g
1928	Mccntot(ig) = Mccntot(ig) +
1929	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
1930	& *(zplev(ig,l) - zplev(ig,l+1)) / g
1931	cccc Column integrated effective ice radius
1932	cccc is weighted by total ice surface area (BETTER than total ice mass)
1933	rave(ig) = rave(ig) +
1934	& tauscaling(ig) *
1935	& zq(ig,l,igcm_ccn_number) *
1936	& (zplev(ig,l) - zplev(ig,l+1)) / g *
1937	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
1938	enddo
1939	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
1940	& /max(pi*rave(ig),1.e-30) ! surface weight
1941	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1942	enddo
1943	else ! of if (scavenging)
1944	rave(:)=0
1945	do ig=1,ngrid
1946	do l=1,nlayer
1947	rave(ig) = rave(ig) +
1948	& zq(ig,l,igcm_h2o_ice) *
1949	& (zplev(ig,l) - zplev(ig,l+1)) / g *
1950	& rice(ig,l) * (1.+nuice_ref)
1951	enddo
1952	rave(ig) = max(rave(ig) /
1953	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
1954	enddo
1955	endif ! of if (scavenging)
1956
1957	!Alternative A. Pottier weighting
1958	rave2(:) = 0.
1959	totrave2(:) = 0.
1960	do ig=1,ngrid
1961	do l=1,nlayer
1962	rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
1963	totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
1964	end do
1965	rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
1966	end do
1967
1968	endif ! of if (water)
1969
1970
1971	endif ! of if (tracer)
1972
1973	#ifndef MESOSCALE
1974	c -----------------------------------------------------------------
1975	c WSTATS: Saving statistics
1976	c -----------------------------------------------------------------
1977	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1978	c which can later be used to make the statistic files of the run:
1979	c "stats") only possible in 3D runs !
1980
1981	IF (callstats) THEN
1982
1983	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1984	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1985	call wstats(ngrid,"co2ice","CO2 ice cover",
1986	& "kg.m-2",2,co2ice)
1987	call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
1988	& 2,tauref)
1989	call wstats(ngrid,"fluxsurf_lw",
1990	& "Thermal IR radiative flux to surface","W.m-2",2,
1991	& fluxsurf_lw)
1992	call wstats(ngrid,"fluxsurf_sw",
1993	& "Solar radiative flux to surface","W.m-2",2,
1994	& fluxsurf_sw_tot)
1995	call wstats(ngrid,"fluxtop_lw",
1996	& "Thermal IR radiative flux to space","W.m-2",2,
1997	& fluxtop_lw)
1998	call wstats(ngrid,"fluxtop_sw",
1999	& "Solar radiative flux to space","W.m-2",2,
2000	& fluxtop_sw_tot)
2001	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2002	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2003	call wstats(ngrid,"v","Meridional (North-South) wind",
2004	& "m.s-1",3,zv)
2005	call wstats(ngrid,"w","Vertical (down-up) wind",
2006	& "m.s-1",3,pw)
2007	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
2008	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
2009	call wstats(ngrid,"q2",
2010	& "Boundary layer eddy kinetic energy",
2011	& "m2.s-2",3,q2)
2012	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
2013	& emis)
2014	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
2015	c & 2,zstress)
2016	c call wstats(ngrid,"sw_htrt","sw heat.rate",
2017	c & "W.m-2",3,zdtsw)
2018	c call wstats(ngrid,"lw_htrt","lw heat.rate",
2019	c & "W.m-2",3,zdtlw)
2020
2021	if (calltherm) then
2022	call wstats(ngrid,"zmax_th","Height of thermals",
2023	& "m",2,zmax_th)
2024	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
2025	& "K.m/s",2,hfmax_th)
2026	call wstats(ngrid,"wstar",
2027	& "Max vertical velocity in thermals",
2028	& "m/s",2,wstar)
2029	endif
2030
2031	if (tracer) then
2032	if (water) then
2033	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2034	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2035	call wstats(ngrid,"vmr_h2ovap",
2036	& "H2O vapor volume mixing ratio","mol/mol",
2037	& 3,vmr)
2038	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2039	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2040	call wstats(ngrid,"vmr_h2oice",
2041	& "H2O ice volume mixing ratio","mol/mol",
2042	& 3,vmr)
2043	! also store vmr_ice*rice for better diagnostics of rice
2044	vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
2045	& rice(1:ngrid,1:nlayer)
2046	call wstats(ngrid,"vmr_h2oice_rice",
2047	& "H2O ice mixing ratio times ice particule size",
2048	& "(mol/mol)*m",
2049	& 3,vmr)
2050	vmr=zqsat(1:ngrid,1:nlayer)
2051	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2052	call wstats(ngrid,"vmr_h2osat",
2053	& "saturation volume mixing ratio","mol/mol",
2054	& 3,vmr)
2055	call wstats(ngrid,"h2o_ice_s",
2056	& "surface h2o_ice","kg/m2",
2057	& 2,qsurf(1,igcm_h2o_ice))
2058	call wstats(ngrid,'albedo',
2059	& 'albedo',
2060	& '',2,albedo(1,1))
2061	call wstats(ngrid,"mtot",
2062	& "total mass of water vapor","kg/m2",
2063	& 2,mtot)
2064	call wstats(ngrid,"icetot",
2065	& "total mass of water ice","kg/m2",
2066	& 2,icetot)
2067	call wstats(ngrid,"reffice",
2068	& "Mean reff","m",
2069	& 2,rave)
2070	call wstats(ngrid,"Nccntot",
2071	& "condensation nuclei","Nbr/m2",
2072	& 2,Nccntot)
2073	call wstats(ngrid,"Mccntot",
2074	& "condensation nuclei mass","kg/m2",
2075	& 2,Mccntot)
2076	call wstats(ngrid,"rice",
2077	& "Ice particle size","m",
2078	& 3,rice)
2079	if (.not.activice) then
2080	call wstats(ngrid,"tauTESap",
2081	& "tau abs 825 cm-1","",
2082	& 2,tauTES)
2083	else
2084	call wstats(ngrid,'tauTES',
2085	& 'tau abs 825 cm-1',
2086	& '',2,taucloudtes)
2087	endif
2088
2089	endif ! of if (water)
2090
2091
2092	if (dustbin.ne.0) then
2093
2094	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
2095
2096	if (doubleq) then
2097	c call wstats(ngrid,'qsurf','qsurf',
2098	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
2099	c call wstats(ngrid,'Nsurf','N particles',
2100	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
2101	c call wstats(ngrid,'dqsdev','ddevil lift',
2102	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
2103	c call wstats(ngrid,'dqssed','sedimentation',
2104	c & 'kg.m-2.s-1',2,zdqssed(1,1))
2105	c call wstats(ngrid,'dqsdif','diffusion',
2106	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
2107	call wstats(ngrid,'dqsdust',
2108	& 'deposited surface dust mass',
2109	& 'kg.m-2.s-1',2,dqdustsurf)
2110	call wstats(ngrid,'dqndust',
2111	& 'deposited surface dust number',
2112	& 'number.m-2.s-1',2,dndustsurf)
2113	call wstats(ngrid,'reffdust','reffdust',
2114	& 'm',3,rdust*ref_r0)
2115	call wstats(ngrid,'dustq','Dust mass mr',
2116	& 'kg/kg',3,qdust)
2117	call wstats(ngrid,'dustN','Dust number',
2118	& 'part/kg',3,ndust)
2119	else
2120	do iq=1,dustbin
2121	write(str2(1:2),'(i2.2)') iq
2122	call wstats(ngrid,'q'//str2,'mix. ratio',
2123	& 'kg/kg',3,zq(1,1,iq))
2124	call wstats(ngrid,'qsurf'//str2,'qsurf',
2125	& 'kg.m-2',2,qsurf(1,iq))
2126	end do
2127	endif ! (doubleq)
2128
2129	if (scavenging) then
2130	call wstats(ngrid,'ccnq','CCN mass mr',
2131	& 'kg/kg',3,qccn)
2132	call wstats(ngrid,'ccnN','CCN number',
2133	& 'part/kg',3,nccn)
2134	endif ! (scavenging)
2135
2136	endif ! (dustbin.ne.0)
2137
2138	if (thermochem .or. photochem) then
2139	do iq=1,nq
2140	if (noms(iq) .ne. "dust_mass" .and.
2141	$ noms(iq) .ne. "dust_number" .and.
2142	$ noms(iq) .ne. "ccn_mass" .and.
2143	$ noms(iq) .ne. "ccn_number" .and.
2144	$ noms(iq) .ne. "ccnco2_mass" .and.
2145	$ noms(iq) .ne. "ccnco2_number") then
2146
2147	! volume mixing ratio
2148
2149	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2150	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
2151
2152	call wstats(ngrid,"vmr_"//trim(noms(iq)),
2153	$ "Volume mixing ratio","mol/mol",3,vmr)
2154	if ((noms(iq).eq."o")
2155	$ .or. (noms(iq).eq."co2")
2156	$ .or. (noms(iq).eq."o3")
2157	$ .or. (noms(iq).eq."ar")
2158	$ .or. (noms(iq).eq."o2")
2159	$ .or. (noms(iq).eq."h2o_vap") ) then
2160	call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
2161	$ "Volume mixing ratio","mol/mol",3,vmr)
2162	end if
2163
2164	! number density (molecule.cm-3)
2165
2166	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2167	& rho(1:ngrid,1:nlayer)n_avog/
2168	& (1000*mmol(iq))
2169
2170	! call wstats(ngrid,"rho_"//trim(noms(iq)),
2171	! $ "Number density","cm-3",3,rhopart)
2172	! call writediagfi(ngrid,"rho_"//trim(noms(iq)),
2173	! $ "Number density","cm-3",3,rhopart)
2174
2175	! vertical column (molecule.cm-2)
2176
2177	do ig = 1,ngrid
2178	colden(ig,iq) = 0.
2179	end do
2180	do l=1,nlayer
2181	do ig=1,ngrid
2182	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
2183	$ *(zplev(ig,l)-zplev(ig,l+1))
2184	$ 6.022e22/(mmol(iq)g)
2185	end do
2186	end do
2187
2188	call wstats(ngrid,"c_"//trim(noms(iq)),
2189	$ "column","mol cm-2",2,colden(1,iq))
2190	call writediagfi(ngrid,"c_"//trim(noms(iq)),
2191	$ "column","mol cm-2",2,colden(1,iq))
2192
2193	! global mass (g)
2194
2195	call planetwide_sumval(colden(:,iq)/6.022e23
2196	$ mmol(iq)1.e4*cell_area(:),mass(iq))
2197
2198	call writediagfi(ngrid,"mass_"//trim(noms(iq)),
2199	$ "global mass","g",0,mass(iq))
2200
2201	end if ! of if (noms(iq) .ne. "dust_mass" ...)
2202	end do ! of do iq=1,nq
2203	end if ! of if (thermochem .or. photochem)
2204
2205	end if ! of if (tracer)
2206
2207	IF(lastcall) THEN
2208	write (,) "Writing stats..."
2209	call mkstats(ierr)
2210	ENDIF
2211
2212	ENDIF !if callstats
2213
2214	c (Store EOF for Mars Climate database software)
2215	IF (calleofdump) THEN
2216	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
2217	ENDIF
2218	#endif
2219	!endif of ifndef MESOSCALE
2220
2221	#ifdef MESOSCALE
2222
2223	!! see comm_wrf.
2224	!! not needed when an array is already in a shared module.
2225	!! --> example : hfmax_th, zmax_th
2226
2227	CALL allocate_comm_wrf(ngrid,nlayer)
2228
2229	!state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s"
2230	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2231	!state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s"
2232	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2233	!state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2"
2234	comm_SWDOWNZ(1:ngrid) = fluxsurf_sw_tot(1:ngrid)
2235	!state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" ""
2236	comm_TAU_DUST(1:ngrid) = tauref(1:ngrid)
2237	!state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m"
2238	comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
2239	!state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2"
2240	IF (igcm_dust_mass .ne. 0) THEN
2241	comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
2242	ELSE
2243	comm_QSURFDUST(1:ngrid) = 0.
2244	ENDIF
2245	!state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic"
2246	comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
2247	!state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2"
2248	comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
2249	!state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
2250	IF (igcm_h2o_ice .ne. 0) THEN
2251	comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
2252	. * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2253	. * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
2254	ELSE
2255	comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
2256	ENDIF
2257	!state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" ""
2258	if (activice) then
2259	comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
2260	else
2261	comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
2262	endif
2263	!state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m"
2264	comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
2265	if (co2clouds) then
2266	comm_RICE(1:ngrid,1:nlayer) = riceco2(1:ngrid,1:nlayer)
2267	endif
2268
2269	!! calculate sensible heat flux in W/m2 for outputs
2270	!! -- the one computed in vdifc is not the real one
2271	!! -- vdifc must have been called
2272	if (.not.callrichsl) then
2273	sensibFlux(1:ngrid) = zflubid(1:ngrid)
2274	. - capcal(1:ngrid)*zdtsdif(1:ngrid)
2275	else
2276	sensibFlux(1:ngrid) =
2277	& (pplay(1:ngrid,1)/(rpt(1:ngrid,1)))cpp
2278	& sqrt(pu(1:ngrid,1)pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
2279	& +(log(1.+0.7wstar(1:ngrid) + 2.3wstar(1:ngrid)2))2)
2280	& zcdh(1:ngrid)(tsurf(1:ngrid)-zh(1:ngrid,1))
2281	endif
2282
2283	#else
2284	#ifndef MESOINI
2285
2286	c ==========================================================
2287	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
2288	c any variable for diagnostic (output with period
2289	c "ecritphy", set in "run.def")
2290	c ==========================================================
2291	c WRITEDIAGFI can ALSO be called from any other subroutines
2292	c for any variables !!
2293	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
2294	c & emis)
2295	c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
2296	c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
2297	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
2298	& tsurf)
2299	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
2300	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
2301	& ,"kg.m-2",2,co2ice)
2302
2303	call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
2304	& "K",2,zt(1,7))
2305	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
2306	& fluxsurf_lw)
2307	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
2308	& fluxsurf_sw_tot)
2309	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
2310	& fluxtop_lw)
2311	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
2312	& fluxtop_sw_tot)
2313	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2314	call WRITEDIAGFI(ngrid,"Time","Time","sols",0,zday)
2315
2316	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
2317	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
2318	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
2319	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
2320	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
2321	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
2322	call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
2323	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
2324	c & zstress)
2325	c call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
2326	c & 'w.m-2',3,zdtsw)
2327	c call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
2328	c & 'w.m-2',3,zdtlw)
2329
2330	if (.not.activice) then
2331	CALL WRITEDIAGFI(ngrid,'tauTESap',
2332	& 'tau abs 825 cm-1',
2333	& '',2,tauTES)
2334	else
2335	CALL WRITEDIAGFI(ngrid,'tauTES',
2336	& 'tau abs 825 cm-1',
2337	& '',2,taucloudtes)
2338	endif
2339
2340	#else
2341	!!! this is to ensure correct initialisation of mesoscale model
2342	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
2343	& tsurf)
2344	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
2345	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
2346	& co2ice)
2347	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2348	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
2349	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
2350	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
2351	& emis)
2352	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
2353	& "K",3,tsoil)
2354	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
2355	& "K",3,inertiedat)
2356	#endif
2357
2358
2359	c ----------------------------------------------------------
2360	c Outputs of the CO2 cycle
2361	c ----------------------------------------------------------
2362
2363	if (tracer.and.(igcm_co2.ne.0)) then
2364	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
2365	! & "kg/kg",2,zq(1,1,igcm_co2))
2366	call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
2367	& "kg/kg",3,zq(1,1,igcm_co2))
2368
2369	! Compute co2 column
2370	co2col(:)=0
2371	do l=1,nlayer
2372	do ig=1,ngrid
2373	co2col(ig)=co2col(ig)+
2374	& zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
2375	enddo
2376	enddo
2377	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
2378	& co2col)
2379	endif ! of if (tracer.and.(igcm_co2.ne.0))
2380
2381
2382	! Output He tracer, if there is one
2383	if (tracer.and.(igcm_he.ne.0)) then
2384	call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
2385	& "kg/kg",3,zq(1,1,igcm_he))
2386	vmr=zq(1:ngrid,1:nlayer,igcm_he)
2387	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
2388	call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
2389	& 'mol/mol',3,vmr)
2390	endif
2391
2392	c ----------------------------------------------------------
2393	c Outputs of the water cycle
2394	c ----------------------------------------------------------
2395	if (tracer) then
2396	if (water) then
2397
2398	#ifdef MESOINI
2399	!!!! waterice = q01, voir readmeteo.F90
2400	call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
2401	& 'kg/kg',3,
2402	& zq(1:ngrid,1:nlayer,igcm_h2o_ice))
2403	!!!! watervapor = q02, voir readmeteo.F90
2404	call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
2405	& 'kg/kg',3,
2406	& zq(1:ngrid,1:nlayer,igcm_h2o_vap))
2407	!!!! surface waterice qsurf02 (voir readmeteo)
2408	call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
2409	& 'kg.m-2',2,
2410	& qsurf(1:ngrid,igcm_h2o_ice))
2411	#endif
2412
2413	CALL WRITEDIAGFI(ngrid,'mtot',
2414	& 'total mass of water vapor',
2415	& 'kg/m2',2,mtot)
2416	CALL WRITEDIAGFI(ngrid,'icetot',
2417	& 'total mass of water ice',
2418	& 'kg/m2',2,icetot)
2419	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2420	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2421	call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
2422	& 'mol/mol',3,vmr)
2423	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2424	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2425	call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
2426	& 'mol/mol',3,vmr)
2427	CALL WRITEDIAGFI(ngrid,'reffice',
2428	& 'Mean reff',
2429	& 'm',2,rave)
2430	!A. Pottier
2431	CALL WRITEDIAGFI(ngrid,'rmoym',
2432	& 'alternative reffice',
2433	& 'm',2,rave2)
2434	call WRITEDIAGFI(ngrid,'saturation',
2435	& 'h2o vap saturation ratio','dimless',3,satu)
2436	if (scavenging) then
2437	CALL WRITEDIAGFI(ngrid,"Nccntot",
2438	& "condensation nuclei","Nbr/m2",
2439	& 2,Nccntot)
2440	CALL WRITEDIAGFI(ngrid,"Mccntot",
2441	& "mass condensation nuclei","kg/m2",
2442	& 2,Mccntot)
2443	endif
2444	call WRITEDIAGFI(ngrid,'rice','Ice particle size',
2445	& 'm',3,rice)
2446	call WRITEDIAGFI(ngrid,'h2o_ice_s',
2447	& 'surface h2o_ice',
2448	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
2449	CALL WRITEDIAGFI(ngrid,'albedo',
2450	& 'albedo',
2451	& '',2,albedo(1,1))
2452	if (tifeedback) then
2453	call WRITEDIAGSOIL(ngrid,"soiltemp",
2454	& "Soil temperature","K",
2455	& 3,tsoil)
2456	call WRITEDIAGSOIL(ngrid,'soilti',
2457	& 'Soil Thermal Inertia',
2458	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
2459	endif
2460	!A. Pottier
2461	if (CLFvarying) then !AP14 nebulosity
2462	call WRITEDIAGFI(ngrid,'totcloudfrac',
2463	& 'Total cloud fraction',
2464	& ' ',2,totcloudfrac)
2465	endif !clf varying
2466
2467	endif !(water)
2468
2469
2470	if (water.and..not.photochem) then
2471	iq=nq
2472	c write(str2(1:2),'(i2.2)') iq
2473	c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
2474	c & 'kg.m-2',2,zdqscloud(1,iq))
2475	c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
2476	c & 'kg/kg',3,zdqchim(1,1,iq))
2477	c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
2478	c & 'kg/kg',3,zdqdif(1,1,iq))
2479	c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
2480	c & 'kg/kg',3,zdqadj(1,1,iq))
2481	c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
2482	c & 'kg/kg',3,zdqc(1,1,iq))
2483	endif !(water.and..not.photochem)
2484	endif
2485
2486	c ----------------------------------------------------------
2487	c Outputs of the dust cycle
2488	c ----------------------------------------------------------
2489
2490	call WRITEDIAGFI(ngrid,'tauref',
2491	& 'Dust ref opt depth','NU',2,tauref)
2492
2493	if (tracer.and.(dustbin.ne.0)) then
2494
2495	call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
2496
2497	#ifndef MESOINI
2498	if (doubleq) then
2499	c call WRITEDIAGFI(ngrid,'qsurf','qsurf',
2500	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
2501	c call WRITEDIAGFI(ngrid,'Nsurf','N particles',
2502	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
2503	c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
2504	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
2505	c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
2506	c & 'kg.m-2.s-1',2,zdqssed(1,1))
2507	c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
2508	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
2509	c call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
2510	c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
2511	c call WRITEDIAGFI(ngrid,'subice','sublimated ice',
2512	c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
2513	call WRITEDIAGFI(ngrid,'dqsdust',
2514	& 'deposited surface dust mass',
2515	& 'kg.m-2.s-1',2,dqdustsurf)
2516	call WRITEDIAGFI(ngrid,'dqndust',
2517	& 'deposited surface dust number',
2518	& 'number.m-2.s-1',2,dndustsurf)
2519	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
2520	& 'm',3,rdust*ref_r0)
2521	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2522	& 'kg/kg',3,qdust)
2523	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2524	& 'part/kg',3,ndust)
2525	call WRITEDIAGFI(ngrid,'dsodust',
2526	& 'dust density scaled opacity',
2527	& 'm2.kg-1',3,dsodust)
2528	c call WRITEDIAGFI(ngrid,"tauscaling",
2529	c & "dust conversion factor"," ",2,tauscaling)
2530	else
2531	do iq=1,dustbin
2532	write(str2(1:2),'(i2.2)') iq
2533	call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
2534	& 'kg/kg',3,zq(1,1,iq))
2535	call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
2536	& 'kg.m-2',2,qsurf(1,iq))
2537	end do
2538	endif ! (doubleq)
2539
2540	if (scavenging) then
2541	call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
2542	& 'kg/kg',3,qccn)
2543	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2544	& 'part/kg',3,nccn)
2545	endif ! (scavenging)
2546
2547	c if (submicron) then
2548	c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
2549	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
2550	c endif ! (submicron)
2551
2552	#else
2553	! !!! to initialize mesoscale we need scaled variables
2554	! !!! because this must correspond to starting point for tracers
2555	! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2556	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
2557	! call WRITEDIAGFI(ngrid,'dustN','Dust number',
2558	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
2559	! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2560	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
2561	! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2562	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
2563	if (freedust) then
2564	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2565	& 'kg/kg',3,qdust)
2566	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2567	& 'part/kg',3,ndust)
2568	call WRITEDIAGFI(ngrid,'ccn','CCN mass mr',
2569	& 'kg/kg',3,qccn)
2570	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2571	& 'part/kg',3,nccn)
2572	else
2573	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2574	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
2575	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2576	& 'part/kg',3,pq(1,1,igcm_dust_number))
2577	call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2578	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
2579	call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2580	& 'part/kg',3,pq(1,1,igcm_ccn_number))
2581	endif
2582	#endif
2583
2584	end if ! (tracer.and.(dustbin.ne.0))
2585
2586
2587	c ----------------------------------------------------------
2588	c Thermospheric outputs
2589	c ----------------------------------------------------------
2590
2591	if(callthermos) then
2592
2593	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
2594	$ 3,zdtnlte)
2595	call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
2596	$ 3,zdteuv)
2597	call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
2598	$ 3,zdtconduc)
2599	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
2600	$ 3,zdtnirco2)
2601
2602	endif !(callthermos)
2603
2604	c ----------------------------------------------------------
2605	c ----------------------------------------------------------
2606	c PBL OUTPUS
2607	c ----------------------------------------------------------
2608	c ----------------------------------------------------------
2609
2610	c ----------------------------------------------------------
2611	c Outputs of thermals
2612	c ----------------------------------------------------------
2613	if (calltherm) then
2614
2615	! call WRITEDIAGFI(ngrid,'dtke',
2616	! & 'tendance tke thermiques','m2/s2',
2617	! & 3,dtke_th)
2618	! call WRITEDIAGFI(ngrid,'d_u_ajs',
2619	! & 'tendance u thermiques','m/s',
2620	! & 3,pdu_th*ptimestep)
2621	! call WRITEDIAGFI(ngrid,'d_v_ajs',
2622	! & 'tendance v thermiques','m/s',
2623	! & 3,pdv_th*ptimestep)
2624	! if (tracer) then
2625	! if (nq .eq. 2) then
2626	! call WRITEDIAGFI(ngrid,'deltaq_th',
2627	! & 'delta q thermiques','kg/kg',
2628	! & 3,ptimestep*pdq_th(:,:,2))
2629	! endif
2630	! endif
2631
2632	call WRITEDIAGFI(ngrid,'zmax_th',
2633	& 'hauteur du thermique','m',
2634	& 2,zmax_th)
2635	call WRITEDIAGFI(ngrid,'hfmax_th',
2636	& 'maximum TH heat flux','K.m/s',
2637	& 2,hfmax_th)
2638	call WRITEDIAGFI(ngrid,'wstar',
2639	& 'maximum TH vertical velocity','m/s',
2640	& 2,wstar)
2641
2642	endif
2643
2644	c ----------------------------------------------------------
2645	c ----------------------------------------------------------
2646	c END OF PBL OUTPUS
2647	c ----------------------------------------------------------
2648	c ----------------------------------------------------------
2649
2650
2651	c ----------------------------------------------------------
2652	c Output in netcdf file "diagsoil.nc" for subterranean
2653	c variables (output every "ecritphy", as for writediagfi)
2654	c ----------------------------------------------------------
2655
2656	! Write soil temperature
2657	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
2658	! & 3,tsoil)
2659	! Write surface temperature
2660	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
2661	! & 2,tsurf)
2662
2663	c ==========================================================
2664	c END OF WRITEDIAGFI
2665	c ==========================================================
2666	#endif
2667	! of ifdef MESOSCALE
2668
2669	ELSE ! if(ngrid.eq.1)
2670
2671	#ifndef MESOSCALE
2672	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
2673	& zls*180./pi,odpref,tauref
2674	c ----------------------------------------------------------------------
2675	c Output in grads file "g1d" (ONLY when using testphys1d)
2676	c (output at every X physical timestep)
2677	c ----------------------------------------------------------------------
2678	c
2679	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
2680	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
2681	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
2682
2683	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
2684	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
2685	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
2686	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
2687
2688	! THERMALS STUFF 1D
2689	if(calltherm) then
2690
2691	call WRITEDIAGFI(ngrid,'lmax_th',
2692	& 'hauteur du thermique','point',
2693	& 0,lmax_th_out)
2694	call WRITEDIAGFI(ngrid,'zmax_th',
2695	& 'hauteur du thermique','m',
2696	& 0,zmax_th)
2697	call WRITEDIAGFI(ngrid,'hfmax_th',
2698	& 'maximum TH heat flux','K.m/s',
2699	& 0,hfmax_th)
2700	call WRITEDIAGFI(ngrid,'wstar',
2701	& 'maximum TH vertical velocity','m/s',
2702	& 0,wstar)
2703
2704	co2col(:)=0.
2705	if (tracer) then
2706	do l=1,nlayer
2707	do ig=1,ngrid
2708	co2col(ig)=co2col(ig)+
2709	& zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
2710	enddo
2711	enddo
2712
2713	end if
2714	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
2715	& ,'kg/m-2',0,co2col)
2716	endif ! of if (calltherm)
2717
2718	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
2719	& ,'m/s',1,pw)
2720	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
2721	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
2722	& tsurf)
2723	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
2724	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
2725
2726	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
2727	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
2728	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
2729	call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
2730	& "K.s-1",1,dtrad)
2731	call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
2732	& 'w.m-2',1,zdtsw)
2733	call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
2734	& 'w.m-2',1,zdtlw)
2735	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
2736	& ,"kg.m-2",0,co2ice)
2737
2738	call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)
2739	do ig=1,ngrid
2740	do l=1,nlayer
2741	satuco2(ig,l) = zq(ig,l,igcm_co2)*
2742	& (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l)
2743
2744	c write(,) "In PHYSIQMOD, pt,zt,time ",pt(ig,l)
2745	c & ,zt(ig,l),ptime
2746	enddo
2747	enddo
2748
2749	c CALL writeg1d(ngrid,nlayer,zt,'temp','K')
2750	c CALL writeg1d(ngrid,nlayer,riceco2,'riceco2','m')
2751	c CALL writeg1d(ngrid,nlayer,satuco2,'satuco2','satu')
2752
2753
2754	c call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
2755	c & satuco2)
2756	c call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
2757	c & ,1,riceco2)
2758	! or output in diagfi.nc (for testphys1d)
2759	call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
2760	call WRITEDIAGFI(ngrid,'temp','Temperature ',
2761	& 'K JA',1,zt)
2762	c call WRITEDIAGFI(ngrid,'temp2','Temperature ',
2763	c & 'K JA2',1,pt)
2764
2765	if(tracer) then
2766	c CALL writeg1d(ngrid,1,tau,'tau','SI')
2767	do iq=1,nq
2768	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
2769	call WRITEDIAGFI(ngrid,trim(noms(iq)),
2770	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
2771	end do
2772	if (doubleq) then
2773	call WRITEDIAGFI(ngrid,'rdust','rdust',
2774	& 'm',1,rdust)
2775	endif
2776	if (water.AND.tifeedback) then
2777	call WRITEDIAGFI(ngrid,"soiltemp",
2778	& "Soil temperature","K",
2779	& 1,tsoil)
2780	call WRITEDIAGFI(ngrid,'soilti',
2781	& 'Soil Thermal Inertia',
2782	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
2783	endif
2784	end if
2785
2786	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
2787	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2788	IF (water) THEN
2789
2790	if (.not.activice) then
2791
2792	tauTES=0
2793	do l=1,nlayer
2794	Qabsice = min(
2795	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
2796	& )
2797	opTES(1,l)= 0.75 * Qabsice *
2798	& zq(1,l,igcm_h2o_ice) *
2799	& (zplev(1,l) - zplev(1,l+1)) / g
2800	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
2801	tauTES=tauTES+ opTES(1,l)
2802	enddo
2803	CALL WRITEDIAGFI(ngrid,'tauTESap',
2804	& 'tau abs 825 cm-1',
2805	& '',0,tauTES)
2806	else
2807
2808	CALL WRITEDIAGFI(ngrid,'tauTES',
2809	& 'tau abs 825 cm-1',
2810	& '',0,taucloudtes)
2811	endif
2812
2813	mtot = 0
2814	icetot = 0
2815	h2otot = qsurf(1,igcm_h2o_ice)
2816
2817	do l=1,nlayer
2818	mtot = mtot + zq(1,l,igcm_h2o_vap)
2819	& * (zplev(1,l) - zplev(1,l+1)) / g
2820	icetot = icetot + zq(1,l,igcm_h2o_ice)
2821	& * (zplev(1,l) - zplev(1,l+1)) / g
2822	end do
2823	h2otot = h2otot+mtot+icetot
2824
2825	CALL WRITEDIAGFI(ngrid,'h2otot',
2826	& 'h2otot',
2827	& 'kg/m2',0,h2otot)
2828	CALL WRITEDIAGFI(ngrid,'mtot',
2829	& 'mtot',
2830	& 'kg/m2',0,mtot)
2831	CALL WRITEDIAGFI(ngrid,'icetot',
2832	& 'icetot',
2833	& 'kg/m2',0,icetot)
2834
2835	if (scavenging) then
2836
2837	rave = 0
2838	do l=1,nlayer
2839	cccc Column integrated effective ice radius
2840	cccc is weighted by total ice surface area (BETTER)
2841	rave = rave + tauscaling(1) *
2842	& zq(1,l,igcm_ccn_number) *
2843	& (zplev(1,l) - zplev(1,l+1)) / g *
2844	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
2845	enddo
2846	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
2847
2848	Nccntot= 0
2849	call watersat(ngrid*nlayer,zt,zplay,zqsat)
2850	do l=1,nlayer
2851	Nccntot = Nccntot +
2852	& zq(1,l,igcm_ccn_number)*tauscaling(1)
2853	& *(zplev(1,l) - zplev(1,l+1)) / g
2854	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
2855	satu(1,l) = (max(satu(1,l),float(1))-1)
2856	! & * zq(1,l,igcm_h2o_vap) *
2857	! & (zplev(1,l) - zplev(1,l+1)) / g
2858	enddo
2859	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
2860	& satu)
2861	CALL WRITEDIAGFI(ngrid,'Nccntot',
2862	& 'Nccntot',
2863	& 'nbr/m2',0,Nccntot)
2864
2865	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
2866	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
2867	call WRITEDIAGFI(ngrid,'zdqsed_dustN'
2868	&,'sedimentation N','Nbr.m-2.s-1',1,
2869	& zdqsed(1,:,igcm_dust_number))
2870
2871	else ! of if (scavenging)
2872
2873	cccc Column integrated effective ice radius
2874	cccc is weighted by total ice mass (LESS GOOD)
2875	rave = 0
2876	do l=1,nlayer
2877	rave = rave + zq(1,l,igcm_h2o_ice)
2878	& * (zplev(1,l) - zplev(1,l+1)) / g
2879	& * rice(1,l) * (1.+nuice_ref)
2880	enddo
2881	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
2882	endif ! of if (scavenging)
2883
2884
2885	CALL WRITEDIAGFI(ngrid,'reffice',
2886	& 'reffice',
2887	& 'm',0,rave)
2888
2889	!Alternative A. Pottier weighting
2890	rave2 = 0.
2891	totrave2 = 0.
2892	do l=1,nlayer
2893	rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l)
2894	totrave2 = totrave2 + zq(1,l,igcm_h2o_ice)
2895	end do
2896	rave2=max(rave2/max(totrave2,1.e-30),1.e-30)
2897	CALL WRITEDIAGFI(ngrid,'rmoym',
2898	& 'reffice',
2899	& 'm',0,rave2)
2900
2901	do iq=1,nq
2902	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
2903	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
2904	end do
2905
2906	call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
2907	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
2908	call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
2909	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
2910	call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
2911	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
2912	& +zdqcloud(1,:,igcm_h2o_vap))
2913
2914	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
2915	& rice)
2916
2917	if (CLFvarying) then
2918	call WRITEDIAGFI(ngrid,'totcloudfrac',
2919	& 'Total cloud fraction',
2920	& ' ',0,totcloudfrac)
2921	endif !clfvarying
2922
2923	ENDIF ! of IF (water)
2924
2925	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2926	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2927
2928
2929	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
2930
2931	do l=2,nlayer-1
2932	tmean=zt(1,l)
2933	if(zt(1,l).ne.zt(1,l-1))
2934	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
2935	zlocal(l)= zlocal(l-1)
2936	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
2937	enddo
2938	zlocal(nlayer)= zlocal(nlayer-1)-
2939	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
2940	& rnew(1,nlayer)*tmean/g
2941	#endif
2942
2943	END IF ! if(ngrid.ne.1)
2944
2945	icount=icount+1
2946
2947	END SUBROUTINE physiq
2948
2949	END MODULE physiq_mod

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

Download in other formats:

Original Format