Context Navigation

physiq_mod.F @ 1649

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

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

Download in other formats:

Original Format