Context Navigation

physiq_mod.F @ 1912

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

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

Download in other formats:

Original Format