Context Navigation

physiq_mod.F @ 1684

Last change on this file since 1684 was 1684, checked in by emillour, 8 years ago
Mars GCM: Add possibility to fix EUV input as E10.7 value and remove previous system (which used parameter solarcondate). The E10.7 value is now set via callphys.def by parameter "fixed_euv_value" which is only used if solvarmod==0. Guidelines for min/ave/max EUV input: fixed_euv_value=80/140/320. EM + FGG
File size: 112.7 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	& igcm_he
24	use comsoil_h, only: inertiedat, ! soil thermal inertia
25	& tsoil, nsoilmx ! number of subsurface layers
26	use geometry_mod, only: longitude, latitude, cell_area
27	use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
28	use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
29	& zthe, z0, albedo_h2o_ice,
30	& frost_albedo_threshold,
31	& tsurf, co2ice, emis,
32	& capcal, fluxgrd, qsurf
33	use comsaison_h, only: dist_sol, declin, mu0, fract
34	use slope_mod, only: theta_sl, psi_sl
35	use conc_mod, only: rnew, cpnew, mmean
36	use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
37	use dimradmars_mod, only: tauscaling, aerosol,
38	& dtrad, fluxrad_sky, fluxrad, albedo,
39	& naerkind
40	use turb_mod, only: q2, wstar, ustar, sensibFlux,
41	& zmax_th, hfmax_th, turb_resolved
42	use planete_h, only: aphelie, periheli, year_day, peri_day,
43	& obliquit
44	USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
45	use param_v4_h, only: nreact,n_avog,
46	& fill_data_thermos, allocate_param_thermos
47	use iono_h, only: allocate_param_iono
48	#ifdef MESOSCALE
49	use comsoil_h, only: mlayer,layer
50	use surfdat_h, only: z0_default
51	use comm_wrf
52	#else
53	use planetwide_mod
54	use phyredem, only: physdem0, physdem1
55	use eofdump_mod, only: eofdump
56	USE vertical_layers_mod, ONLY: ap,bp,aps,bps
57	#endif
58
59	IMPLICIT NONE
60	c=======================================================================
61	c
62	c subject:
63	c --------
64	c
65	c Organisation of the physical parametrisations of the LMD
66	c martian atmospheric general circulation model.
67	c
68	c The GCM can be run without or with tracer transport
69	c depending on the value of Logical "tracer" in file "callphys.def"
70	c Tracers may be water vapor, ice OR chemical species OR dust particles
71	c
72	c SEE comments in initracer.F about numbering of tracer species...
73	c
74	c It includes:
75	c
76	c 1. Initialization:
77	c 1.1 First call initializations
78	c 1.2 Initialization for every call to physiq
79	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
80	c 2. Compute radiative transfer tendencies
81	c (longwave and shortwave) for CO2 and aerosols.
82	c 3. Gravity wave and subgrid scale topography drag :
83	c 4. Vertical diffusion (turbulent mixing):
84	c 5. Convective adjustment
85	c 6. Condensation and sublimation of carbon dioxide.
86	c 7. TRACERS :
87	c 7a. water, water ice, co2 ice (clouds)
88	c 7b. call for photochemistry when tracers are chemical species
89	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
90	c 7d. updates (CO2 pressure variations, surface budget)
91	c 8. Contribution to tendencies due to thermosphere
92	c 9. Surface and sub-surface temperature calculations
93	c 10. Write outputs :
94	c - "startfi", "histfi" (if it's time)
95	c - Saving statistics (if "callstats = .true.")
96	c - Dumping eof (if "calleofdump = .true.")
97	c - Output any needed variables in "diagfi"
98	c 11. Diagnostic: mass conservation of tracers
99	c
100	c author:
101	c -------
102	c Frederic Hourdin 15/10/93
103	c Francois Forget 1994
104	c Christophe Hourdin 02/1997
105	c Subroutine completly rewritten by F.Forget (01/2000)
106	c Introduction of the photochemical module: S. Lebonnois (11/2002)
107	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
108	c Water ice clouds: Franck Montmessin (update 06/2003)
109	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
110	c Nb: See callradite.F for more information.
111	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
112	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
113	c
114	c 10/16 J. Audouard: modifications for CO2 clouds scheme
115
116	c arguments:
117	c ----------
118	c
119	c input:
120	c ------
121	c ecri period (in dynamical timestep) to write output
122	c ngrid Size of the horizontal grid.
123	c All internal loops are performed on that grid.
124	c nlayer Number of vertical layers.
125	c nq Number of advected fields
126	c firstcall True at the first call
127	c lastcall True at the last call
128	c pday Number of days counted from the North. Spring
129	c equinoxe.
130	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
131	c ptimestep timestep (s)
132	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
133	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
134	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
135	c pu(ngrid,nlayer) u component of the wind (ms-1)
136	c pv(ngrid,nlayer) v component of the wind (ms-1)
137	c pt(ngrid,nlayer) Temperature (K)
138	c pq(ngrid,nlayer,nq) Advected fields
139	c pudyn(ngrid,nlayer) \|
140	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
141	c ptdyn(ngrid,nlayer) \| corresponding variables
142	c pqdyn(ngrid,nlayer,nq) \|
143	c flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
144	c
145	c output:
146	c -------
147	c
148	c pdu(ngrid,nlayer) \|
149	c pdv(ngrid,nlayer) \| Temporal derivative of the corresponding
150	c pdt(ngrid,nlayer) \| variables due to physical processes.
151	c pdq(ngrid,nlayer,nq) \|
152	c pdpsrf(ngrid) \|
153
154	c
155	c=======================================================================
156	c
157	c 0. Declarations :
158	c ------------------
159
160	#include "callkeys.h"
161	#include "comg1d.h"
162	#include "nlteparams.h"
163	#include "chimiedata.h"
164	#include "netcdf.inc"
165
166	c Arguments :
167	c -----------
168
169	c inputs:
170	c -------
171	INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
172	INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
173	INTEGER,INTENT(in) :: nq ! number of tracers
174	LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
175	LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
176	REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
177	REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
178	REAL,INTENT(in) :: ptimestep ! physics timestep (s)
179	REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
180	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
181	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
182	REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
183	REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
184	REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
185	REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
186	REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
187	! at lower boundary of layer
188
189	c outputs:
190	c --------
191	c physical tendencies
192	REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
193	REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
194	REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
195	REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
196	REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
197
198
199	c Local saved variables:
200	c ----------------------
201	INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
202	INTEGER,SAVE :: icount ! counter of calls to physiq during the run.
203	#ifdef DUSTSTORM
204	REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
205	#endif
206	c Variables used by the water ice microphysical scheme:
207	REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m)
208	REAL nuice(ngrid,nlayer) ! Estimated effective variance
209	! of the size distribution
210	real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
211	real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3)
212	REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
213	! thermal inertia (J.s-1/2.m-2.K-1)
214	! (used only when tifeedback=.true.)
215	c Variables used by the CO2 clouds microphysical scheme:
216	REAL riceco2(ngrid,nlayer) ! co2 ice geometric mean radius (m)
217	real rsedcloudco2(ngrid,nlayer) !CO2 Cloud sedimentation radius
218	real rhocloudco2(ngrid,nlayer) !co2 Cloud density (kg.m-3)
219	real zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
220	c Variables used by the photochemistry
221	logical :: asis ! true : adaptative semi-implicit symmetric (asis) chemical solver
222	! false : euler backward chemical solver
223	REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
224	REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry)
225	c Variables used by the slope model
226	REAL sl_ls, sl_lct, sl_lat
227	REAL sl_tau, sl_alb, sl_the, sl_psi
228	REAL sl_fl0, sl_flu
229	REAL sl_ra, sl_di0
230	REAL sky
231
232	REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
233
234	c Local variables :
235	c -----------------
236
237	REAL CBRT
238	EXTERNAL CBRT
239
240	! CHARACTER*80 fichier
241	INTEGER l,ig,ierr,igout,iq,tapphys
242
243	REAL fluxsurf_lw(ngrid) !incident LW (IR) surface flux (W.m-2)
244	REAL fluxsurf_sw(ngrid,2) !incident SW (solar) surface flux (W.m-2)
245	REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2)
246	REAL fluxtop_sw(ngrid,2) !Outgoing SW (solar) flux to space (W.m-2)
247	REAL tauref(ngrid) ! Reference column optical depth at odpref
248	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
249	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
250	! AS: TBD: this one should be in a module !
251	REAL dsodust(ngrid,nlayer) ! density-scaled opacity (in infrared)
252	REAL zls ! solar longitude (rad)
253	REAL zday ! date (time since Ls=0, in martian days)
254	REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers
255	REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
256	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
257
258	c Tendancies due to various processes:
259	REAL dqsurf(ngrid,nq)
260	REAL zdtlw(ngrid,nlayer) ! (K/s)
261	REAL zdtsw(ngrid,nlayer) ! (K/s)
262	! REAL cldtlw(ngrid,nlayer) ! (K/s) LW heating rate for clear area
263	! REAL cldtsw(ngrid,nlayer) ! (K/s) SW heating rate for clear area
264	REAL zdtnirco2(ngrid,nlayer) ! (K/s)
265	REAL zdtnlte(ngrid,nlayer) ! (K/s)
266	REAL zdtsurf(ngrid) ! (K/s)
267	REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
268	REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
269	REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid) ! (K/s)
270	REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
271	REAL zdhadj(ngrid,nlayer) ! (K/s)
272	REAL zdtgw(ngrid,nlayer) ! (K/s)
273	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
274	REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
275	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
276
277	REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
278	REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
279	REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
280	REAL zdqadj(ngrid,nlayer,nq)
281	REAL zdqc(ngrid,nlayer,nq)
282	REAL zdqcloud(ngrid,nlayer,nq),zdqcloudco2(ngrid,nlayer,nq)
283	REAL zdqscloud(ngrid,nq)
284	REAL zdqchim(ngrid,nlayer,nq)
285	REAL zdqschim(ngrid,nq)
286
287	REAL zdteuv(ngrid,nlayer) ! (K/s)
288	REAL zdtconduc(ngrid,nlayer) ! (K/s)
289	REAL zdumolvis(ngrid,nlayer)
290	REAL zdvmolvis(ngrid,nlayer)
291	real zdqmoldiff(ngrid,nlayer,nq)
292
293	c Local variable for local intermediate calcul:
294	REAL zflubid(ngrid)
295	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
296	REAL zdum1(ngrid,nlayer)
297	REAL zdum2(ngrid,nlayer)
298	REAL ztim1,ztim2,ztim3, z1,z2
299	REAL ztime_fin
300	REAL zdh(ngrid,nlayer)
301	REAL zh(ngrid,nlayer) ! potential temperature (K)
302	REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
303	INTEGER length
304	PARAMETER (length=100)
305
306	c local variables only used for diagnostic (output in file "diagfi" or "stats")
307	c -----------------------------------------------------------------------------
308	REAL ps(ngrid), zt(ngrid,nlayer)
309	REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
310	REAL zq(ngrid,nlayer,nq)
311	REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
312	character*2 str2
313	! character*5 str5
314	real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
315	real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
316	integer igmin, lmin
317	logical tdiag
318
319	real co2col(ngrid) ! CO2 column
320	! pplev and pplay are dynamical inputs and must not be modified in the physics.
321	! instead, use zplay and zplev :
322	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
323	! REAL zstress(ngrid),cd
324	real tmean, zlocal(nlayer)
325	real rho(ngrid,nlayer) ! density
326	real vmr(ngrid,nlayer) ! volume mixing ratio
327	real rhopart(ngrid,nlayer) ! number density of a given species
328	real colden(ngrid,nq) ! vertical column of tracers
329	real mass(nq) ! global mass of tracers (g)
330	REAL mtot(ngrid) ! Total mass of water vapor (kg/m2)
331	REAL icetot(ngrid) ! Total mass of water ice (kg/m2)
332	REAL mtotco2(ngrid) ! Total mass of co2 vapor (kg/m2)
333	REAL icetotco2(ngrid) ! Total mass of co2 ice (kg/m2)
334	REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2)
335	REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2)
336	REAL rave(ngrid) ! Mean water ice effective radius (m)
337	REAL raveco2(ngrid) ! Mean co2 ice effective radius (m)
338	REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1
339	REAL tauTES(ngrid) ! column optical depth at 825 cm-1
340	REAL Qabsice ! Water ice absorption coefficient
341	REAL taucloudtes(ngrid) ! Cloud opacity at infrared
342	! reference wavelength using
343	! Qabs instead of Qext
344	! (direct comparison with TES)
345
346	REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
347	REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
348	REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
349	REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
350	REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg)
351	REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg)
352
353	c Test 1d/3d scavenging
354	real h2otot(ngrid)
355	REAL satu(ngrid,nlayer) ! satu ratio for output
356	REAL zqsat(ngrid,nlayer) ! saturation
357	REAL satuco2(ngrid,nlayer) ! co2 satu ratio for output
358	REAL zqsatco2(ngrid,nlayer) ! saturation co2
359	REAL,SAVE :: time_phys
360
361	! Added for new NLTE scheme
362
363	real co2vmr_gcm(ngrid,nlayer)
364	real n2vmr_gcm(ngrid,nlayer)
365	real ovmr_gcm(ngrid,nlayer)
366	real covmr_gcm(ngrid,nlayer)
367	integer ierr_nlte
368	real*8 varerr
369
370	c Variables for PBL
371	REAL zz1(ngrid)
372	REAL lmax_th_out(ngrid)
373	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
374	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
375	INTEGER lmax_th(ngrid),dimout,n_out,n
376	CHARACTER(50) zstring
377	REAL dtke_th(ngrid,nlayer+1)
378	REAL zcdv(ngrid), zcdh(ngrid)
379	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
380	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
381	REAL u_out1(ngrid)
382	REAL T_out1(ngrid)
383	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
384	REAL tstar(ngrid) ! friction velocity and friction potential temp
385	REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
386	! REAL zu2(ngrid)
387
388	c=======================================================================
389
390	c 1. Initialisation:
391	c -----------------
392
393	c 1.1 Initialisation only at first call
394	c ---------------------------------------
395	IF (firstcall) THEN
396
397	c variables set to 0
398	c ~~~~~~~~~~~~~~~~~~
399	aerosol(:,:,:)=0
400	dtrad(:,:)=0
401
402	#ifndef MESOSCALE
403	fluxrad(:)=0
404	wstar(:)=0.
405	#endif
406
407	c read startfi
408	c ~~~~~~~~~~~~
409	#ifndef MESOSCALE
410	! GCM. Read netcdf initial physical parameters.
411	CALL phyetat0 ("startfi.nc",0,0,
412	& nsoilmx,ngrid,nlayer,nq,
413	& day_ini,time_phys,
414	& tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
415
416	if (pday.ne.day_ini) then
417	write(,) "PHYSIQ: ERROR: bad synchronization between ",
418	& "physics and dynamics"
419	write(,) "dynamics day: ",pday
420	write(,) "physics day: ",day_ini
421	stop
422	endif
423
424	write (,) 'In physiq day_ini =', day_ini
425
426	#else
427	! MESOSCALE. Supposedly everything is already set in modules.
428	! So we just check. And we fill day_ini
429	print*,"check: --- in physiq.F"
430	print*,"check: rad,cpp,g,r,rcp,daysec"
431	print*,rad,cpp,g,r,rcp,daysec
432	PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
433	PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
434	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
435	PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
436	PRINT*,'check: tracernames ', noms
437	PRINT*,'check: emis ',emis(1),emis(ngrid)
438	PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
439	PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
440	PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid)
441	!!!
442	day_ini = pday
443	#endif
444
445	c initialize tracers
446	c ~~~~~~~~~~~~~~~~~~
447	IF (tracer) THEN
448	CALL initracer(ngrid,nq,qsurf)
449	ENDIF ! end tracer
450
451	c Initialize albedo and orbital calculation
452	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453	CALL surfini(ngrid,co2ice,qsurf,albedo)
454	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
455
456	c initialize soil
457	c ~~~~~~~~~~~~~~~
458	IF (callsoil) THEN
459	c Thermal inertia feedback:
460	IF (tifeedback) THEN
461	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
462	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
463	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
464	ELSE
465	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
466	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
467	ENDIF ! of IF (tifeedback)
468	ELSE
469	PRINT*,
470	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
471	DO ig=1,ngrid
472	capcal(ig)=1.e5
473	fluxgrd(ig)=0.
474	ENDDO
475	ENDIF
476	icount=1
477
478	#ifndef MESOSCALE
479	c Initialize thermospheric parameters
480	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
481
482	if (callthermos) then
483	call fill_data_thermos
484	call allocate_param_thermos(nlayer)
485	call allocate_param_iono(nlayer,nreact)
486	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
2339	! Output He tracer, if there is one
2340	if (tracer.and.(igcm_he.ne.0)) then
2341	call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
2342	& "kg/kg",3,zq(1,1,igcm_he))
2343	vmr=zq(1:ngrid,1:nlayer,igcm_he)
2344	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
2345	call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
2346	& 'mol/mol',3,vmr)
2347	endif
2348
2349	c ----------------------------------------------------------
2350	c Outputs of the water cycle
2351	c ----------------------------------------------------------
2352	if (tracer) then
2353	if (water) then
2354
2355	#ifdef MESOINI
2356	!!!! waterice = q01, voir readmeteo.F90
2357	call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
2358	& 'kg/kg',3,
2359	& zq(1:ngrid,1:nlayer,igcm_h2o_ice))
2360	!!!! watervapor = q02, voir readmeteo.F90
2361	call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
2362	& 'kg/kg',3,
2363	& zq(1:ngrid,1:nlayer,igcm_h2o_vap))
2364	!!!! surface waterice qsurf02 (voir readmeteo)
2365	call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
2366	& 'kg.m-2',2,
2367	& qsurf(1:ngrid,igcm_h2o_ice))
2368	#endif
2369
2370	CALL WRITEDIAGFI(ngrid,'mtot',
2371	& 'total mass of water vapor',
2372	& 'kg/m2',2,mtot)
2373	CALL WRITEDIAGFI(ngrid,'icetot',
2374	& 'total mass of water ice',
2375	& 'kg/m2',2,icetot)
2376	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2377	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2378	call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
2379	& 'mol/mol',3,vmr)
2380	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2381	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2382	call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
2383	& 'mol/mol',3,vmr)
2384	CALL WRITEDIAGFI(ngrid,'reffice',
2385	& 'Mean reff',
2386	& 'm',2,rave)
2387	call WRITEDIAGFI(ngrid,'saturation',
2388	& 'h2o vap saturation ratio','dimless',3,satu)
2389	if (scavenging) then
2390	CALL WRITEDIAGFI(ngrid,"Nccntot",
2391	& "condensation nuclei","Nbr/m2",
2392	& 2,Nccntot)
2393	CALL WRITEDIAGFI(ngrid,"Mccntot",
2394	& "mass condensation nuclei","kg/m2",
2395	& 2,Mccntot)
2396	endif
2397	call WRITEDIAGFI(ngrid,'rice','Ice particle size',
2398	& 'm',3,rice)
2399	call WRITEDIAGFI(ngrid,'h2o_ice_s',
2400	& 'surface h2o_ice',
2401	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
2402	CALL WRITEDIAGFI(ngrid,'albedo',
2403	& 'albedo',
2404	& '',2,albedo(1,1))
2405	if (tifeedback) then
2406	call WRITEDIAGSOIL(ngrid,"soiltemp",
2407	& "Soil temperature","K",
2408	& 3,tsoil)
2409	call WRITEDIAGSOIL(ngrid,'soilti',
2410	& 'Soil Thermal Inertia',
2411	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
2412	endif
2413	endif !(water)
2414
2415
2416	if (water.and..not.photochem) then
2417	iq=nq
2418	c write(str2(1:2),'(i2.2)') iq
2419	c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
2420	c & 'kg.m-2',2,zdqscloud(1,iq))
2421	c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
2422	c & 'kg/kg',3,zdqchim(1,1,iq))
2423	c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
2424	c & 'kg/kg',3,zdqdif(1,1,iq))
2425	c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
2426	c & 'kg/kg',3,zdqadj(1,1,iq))
2427	c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
2428	c & 'kg/kg',3,zdqc(1,1,iq))
2429	endif !(water.and..not.photochem)
2430	endif
2431
2432	c ----------------------------------------------------------
2433	c Outputs of the dust cycle
2434	c ----------------------------------------------------------
2435
2436	call WRITEDIAGFI(ngrid,'tauref',
2437	& 'Dust ref opt depth','NU',2,tauref)
2438
2439	if (tracer.and.(dustbin.ne.0)) then
2440
2441	call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
2442
2443	#ifndef MESOINI
2444	if (doubleq) then
2445	c call WRITEDIAGFI(ngrid,'qsurf','qsurf',
2446	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
2447	c call WRITEDIAGFI(ngrid,'Nsurf','N particles',
2448	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
2449	c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
2450	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
2451	c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
2452	c & 'kg.m-2.s-1',2,zdqssed(1,1))
2453	c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
2454	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
2455	c call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
2456	c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
2457	c call WRITEDIAGFI(ngrid,'subice','sublimated ice',
2458	c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
2459	call WRITEDIAGFI(ngrid,'dqsdust',
2460	& 'deposited surface dust mass',
2461	& 'kg.m-2.s-1',2,dqdustsurf)
2462	call WRITEDIAGFI(ngrid,'dqndust',
2463	& 'deposited surface dust number',
2464	& 'number.m-2.s-1',2,dndustsurf)
2465	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
2466	& 'm',3,rdust*ref_r0)
2467	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2468	& 'kg/kg',3,qdust)
2469	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2470	& 'part/kg',3,ndust)
2471	call WRITEDIAGFI(ngrid,'dsodust',
2472	& 'dust density scaled opacity',
2473	& 'm2.kg-1',3,dsodust)
2474	c call WRITEDIAGFI(ngrid,"tauscaling",
2475	c & "dust conversion factor"," ",2,tauscaling)
2476	else
2477	do iq=1,dustbin
2478	write(str2(1:2),'(i2.2)') iq
2479	call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
2480	& 'kg/kg',3,zq(1,1,iq))
2481	call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
2482	& 'kg.m-2',2,qsurf(1,iq))
2483	end do
2484	endif ! (doubleq)
2485
2486	if (scavenging) then
2487	call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
2488	& 'kg/kg',3,qccn)
2489	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2490	& 'part/kg',3,nccn)
2491	endif ! (scavenging)
2492
2493	c if (submicron) then
2494	c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
2495	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
2496	c endif ! (submicron)
2497
2498	#else
2499	! !!! to initialize mesoscale we need scaled variables
2500	! !!! because this must correspond to starting point for tracers
2501	! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2502	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
2503	! call WRITEDIAGFI(ngrid,'dustN','Dust number',
2504	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
2505	! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2506	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
2507	! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2508	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
2509	if (freedust) then
2510	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2511	& 'kg/kg',3,qdust)
2512	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2513	& 'part/kg',3,ndust)
2514	call WRITEDIAGFI(ngrid,'ccn','CCN mass mr',
2515	& 'kg/kg',3,qccn)
2516	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2517	& 'part/kg',3,nccn)
2518	else
2519	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2520	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
2521	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2522	& 'part/kg',3,pq(1,1,igcm_dust_number))
2523	call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2524	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
2525	call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2526	& 'part/kg',3,pq(1,1,igcm_ccn_number))
2527	endif
2528	#endif
2529
2530	end if ! (tracer.and.(dustbin.ne.0))
2531
2532
2533	c ----------------------------------------------------------
2534	c Thermospheric outputs
2535	c ----------------------------------------------------------
2536
2537	if(callthermos) then
2538
2539	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
2540	$ 3,zdtnlte)
2541	call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
2542	$ 3,zdteuv)
2543	call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
2544	$ 3,zdtconduc)
2545	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
2546	$ 3,zdtnirco2)
2547
2548	endif !(callthermos)
2549
2550	c ----------------------------------------------------------
2551	c ----------------------------------------------------------
2552	c PBL OUTPUS
2553	c ----------------------------------------------------------
2554	c ----------------------------------------------------------
2555
2556	c ----------------------------------------------------------
2557	c Outputs of thermals
2558	c ----------------------------------------------------------
2559	if (calltherm) then
2560
2561	! call WRITEDIAGFI(ngrid,'dtke',
2562	! & 'tendance tke thermiques','m2/s2',
2563	! & 3,dtke_th)
2564	! call WRITEDIAGFI(ngrid,'d_u_ajs',
2565	! & 'tendance u thermiques','m/s',
2566	! & 3,pdu_th*ptimestep)
2567	! call WRITEDIAGFI(ngrid,'d_v_ajs',
2568	! & 'tendance v thermiques','m/s',
2569	! & 3,pdv_th*ptimestep)
2570	! if (tracer) then
2571	! if (nq .eq. 2) then
2572	! call WRITEDIAGFI(ngrid,'deltaq_th',
2573	! & 'delta q thermiques','kg/kg',
2574	! & 3,ptimestep*pdq_th(:,:,2))
2575	! endif
2576	! endif
2577
2578	call WRITEDIAGFI(ngrid,'zmax_th',
2579	& 'hauteur du thermique','m',
2580	& 2,zmax_th)
2581	call WRITEDIAGFI(ngrid,'hfmax_th',
2582	& 'maximum TH heat flux','K.m/s',
2583	& 2,hfmax_th)
2584	call WRITEDIAGFI(ngrid,'wstar',
2585	& 'maximum TH vertical velocity','m/s',
2586	& 2,wstar)
2587
2588	endif
2589
2590	c ----------------------------------------------------------
2591	c ----------------------------------------------------------
2592	c END OF PBL OUTPUS
2593	c ----------------------------------------------------------
2594	c ----------------------------------------------------------
2595
2596
2597	c ----------------------------------------------------------
2598	c Output in netcdf file "diagsoil.nc" for subterranean
2599	c variables (output every "ecritphy", as for writediagfi)
2600	c ----------------------------------------------------------
2601
2602	! Write soil temperature
2603	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
2604	! & 3,tsoil)
2605	! Write surface temperature
2606	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
2607	! & 2,tsurf)
2608
2609	c ==========================================================
2610	c END OF WRITEDIAGFI
2611	c ==========================================================
2612	#endif
2613	! of ifdef MESOSCALE
2614
2615	ELSE ! if(ngrid.eq.1)
2616
2617	#ifndef MESOSCALE
2618	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
2619	& zls*180./pi,odpref,tauref
2620	c ----------------------------------------------------------------------
2621	c Output in grads file "g1d" (ONLY when using testphys1d)
2622	c (output at every X physical timestep)
2623	c ----------------------------------------------------------------------
2624	c
2625	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
2626	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
2627	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
2628
2629	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
2630	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
2631	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
2632	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
2633
2634	! THERMALS STUFF 1D
2635	if(calltherm) then
2636
2637	call WRITEDIAGFI(ngrid,'lmax_th',
2638	& 'hauteur du thermique','point',
2639	& 0,lmax_th_out)
2640	call WRITEDIAGFI(ngrid,'zmax_th',
2641	& 'hauteur du thermique','m',
2642	& 0,zmax_th)
2643	call WRITEDIAGFI(ngrid,'hfmax_th',
2644	& 'maximum TH heat flux','K.m/s',
2645	& 0,hfmax_th)
2646	call WRITEDIAGFI(ngrid,'wstar',
2647	& 'maximum TH vertical velocity','m/s',
2648	& 0,wstar)
2649
2650	co2col(:)=0.
2651	if (tracer) then
2652	do l=1,nlayer
2653	do ig=1,ngrid
2654	co2col(ig)=co2col(ig)+
2655	& zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
2656	enddo
2657	enddo
2658
2659	end if
2660	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
2661	& ,'kg/m-2',0,co2col)
2662	endif ! of if (calltherm)
2663
2664	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
2665	& ,'m/s',1,pw)
2666	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
2667	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
2668	& tsurf)
2669	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
2670	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
2671
2672	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
2673	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
2674	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
2675	call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
2676	& "K.s-1",1,dtrad)
2677	call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
2678	& 'w.m-2',1,zdtsw)
2679	call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
2680	& 'w.m-2',1,zdtlw)
2681	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
2682	& ,"kg.m-2",0,co2ice)
2683
2684	call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)
2685	do ig=1,ngrid
2686	do l=1,nlayer
2687	satuco2(ig,l) = zq(ig,l,igcm_co2)*
2688	& (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l)
2689
2690	c write(,) "In PHYSIQMOD, pt,zt,time ",pt(ig,l)
2691	c & ,zt(ig,l),ptime
2692	enddo
2693	enddo
2694
2695	c CALL writeg1d(ngrid,nlayer,zt,'temp','K')
2696	c CALL writeg1d(ngrid,nlayer,riceco2,'riceco2','m')
2697	c CALL writeg1d(ngrid,nlayer,satuco2,'satuco2','satu')
2698
2699
2700	c call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
2701	c & satuco2)
2702	c call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
2703	c & ,1,riceco2)
2704	! or output in diagfi.nc (for testphys1d)
2705	call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
2706	call WRITEDIAGFI(ngrid,'temp','Temperature ',
2707	& 'K JA',1,zt)
2708	c call WRITEDIAGFI(ngrid,'temp2','Temperature ',
2709	c & 'K JA2',1,pt)
2710
2711	if(tracer) then
2712	c CALL writeg1d(ngrid,1,tau,'tau','SI')
2713	do iq=1,nq
2714	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
2715	call WRITEDIAGFI(ngrid,trim(noms(iq)),
2716	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
2717	end do
2718	if (doubleq) then
2719	call WRITEDIAGFI(ngrid,'rdust','rdust',
2720	& 'm',1,rdust)
2721	endif
2722	if (water.AND.tifeedback) then
2723	call WRITEDIAGFI(ngrid,"soiltemp",
2724	& "Soil temperature","K",
2725	& 1,tsoil)
2726	call WRITEDIAGFI(ngrid,'soilti',
2727	& 'Soil Thermal Inertia',
2728	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
2729	endif
2730	end if
2731
2732	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
2733	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2734	IF (water) THEN
2735
2736	if (.not.activice) then
2737
2738	tauTES=0
2739	do l=1,nlayer
2740	Qabsice = min(
2741	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
2742	& )
2743	opTES(1,l)= 0.75 * Qabsice *
2744	& zq(1,l,igcm_h2o_ice) *
2745	& (zplev(1,l) - zplev(1,l+1)) / g
2746	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
2747	tauTES=tauTES+ opTES(1,l)
2748	enddo
2749	CALL WRITEDIAGFI(ngrid,'tauTESap',
2750	& 'tau abs 825 cm-1',
2751	& '',0,tauTES)
2752	else
2753
2754	CALL WRITEDIAGFI(ngrid,'tauTES',
2755	& 'tau abs 825 cm-1',
2756	& '',0,taucloudtes)
2757	endif
2758
2759	mtot = 0
2760	icetot = 0
2761	h2otot = qsurf(1,igcm_h2o_ice)
2762
2763	do l=1,nlayer
2764	mtot = mtot + zq(1,l,igcm_h2o_vap)
2765	& * (zplev(1,l) - zplev(1,l+1)) / g
2766	icetot = icetot + zq(1,l,igcm_h2o_ice)
2767	& * (zplev(1,l) - zplev(1,l+1)) / g
2768	end do
2769	h2otot = h2otot+mtot+icetot
2770
2771	CALL WRITEDIAGFI(ngrid,'h2otot',
2772	& 'h2otot',
2773	& 'kg/m2',0,h2otot)
2774	CALL WRITEDIAGFI(ngrid,'mtot',
2775	& 'mtot',
2776	& 'kg/m2',0,mtot)
2777	CALL WRITEDIAGFI(ngrid,'icetot',
2778	& 'icetot',
2779	& 'kg/m2',0,icetot)
2780
2781	if (scavenging) then
2782
2783	rave = 0
2784	do l=1,nlayer
2785	cccc Column integrated effective ice radius
2786	cccc is weighted by total ice surface area (BETTER)
2787	rave = rave + tauscaling(1) *
2788	& zq(1,l,igcm_ccn_number) *
2789	& (zplev(1,l) - zplev(1,l+1)) / g *
2790	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
2791	enddo
2792	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
2793
2794	Nccntot= 0
2795	call watersat(ngrid*nlayer,zt,zplay,zqsat)
2796	do l=1,nlayer
2797	Nccntot = Nccntot +
2798	& zq(1,l,igcm_ccn_number)*tauscaling(1)
2799	& *(zplev(1,l) - zplev(1,l+1)) / g
2800	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
2801	satu(1,l) = (max(satu(1,l),float(1))-1)
2802	! & * zq(1,l,igcm_h2o_vap) *
2803	! & (zplev(1,l) - zplev(1,l+1)) / g
2804	enddo
2805	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
2806	& satu)
2807	CALL WRITEDIAGFI(ngrid,'Nccntot',
2808	& 'Nccntot',
2809	& 'nbr/m2',0,Nccntot)
2810
2811	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
2812	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
2813	call WRITEDIAGFI(ngrid,'zdqsed_dustN'
2814	&,'sedimentation N','Nbr.m-2.s-1',1,
2815	& zdqsed(1,:,igcm_dust_number))
2816
2817	else ! of if (scavenging)
2818
2819	cccc Column integrated effective ice radius
2820	cccc is weighted by total ice mass (LESS GOOD)
2821	rave = 0
2822	do l=1,nlayer
2823	rave = rave + zq(1,l,igcm_h2o_ice)
2824	& * (zplev(1,l) - zplev(1,l+1)) / g
2825	& * rice(1,l) * (1.+nuice_ref)
2826	enddo
2827	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
2828	endif ! of if (scavenging)
2829
2830	if (co2clouds) then
2831	if (scavenging) then
2832	Nccntot(:)= 0
2833	Mccntot(:)= 0
2834	raveco2(:)=0
2835	icetotco2(:)=0
2836	do ig=1,ngrid
2837	do l=1,nlayer
2838	icetotco2(ig) = icetotco2(ig) +
2839	& zq(ig,l,igcm_co2_ice) *
2840	& (zplev(ig,l) - zplev(ig,l+1)) / g
2841	Nccntot(ig) = Nccntot(ig) +
2842	& zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)
2843	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2844	Mccntot(ig) = Mccntot(ig) +
2845	& zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig)
2846	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2847	cccc Column integrated effective ice radius
2848	cccc is weighted by total ice surface area (BETTER than total ice mass)
2849	raveco2(ig) = raveco2(ig) +
2850	& tauscaling(ig) *
2851	& zq(ig,l,igcm_ccnco2_number) *
2852	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2853	& riceco2(ig,l) * riceco2(ig,l)*
2854	& (1.+nuiceco2_ref)
2855	enddo
2856	raveco2(ig)=(icetotco2(ig)/rho_ice_co2+Mccntot(ig)/
2857	& rho_dust)*0.75
2858	& /max(pi*raveco2(ig),1.e-30) ! surface weight
2859	if (icetotco2(ig)*1e3.lt.0.01) raveco2(ig)=0.
2860	enddo
2861	else ! of if (scavenging)
2862	raveco2(:)=0
2863	do ig=1,ngrid
2864	do l=1,nlayer
2865	raveco2(ig) = raveco2(ig) +
2866	& zq(ig,l,igcm_co2_ice) *
2867	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2868	& riceco2(ig,l) * (1.+nuiceco2_ref)
2869	enddo
2870	raveco2(ig) = max(raveco2(ig) /
2871	& max(icetotco2(ig),1.e-30),1.e-30) ! mass weight
2872	enddo
2873	endif ! of if (scavenging)
2874
2875	call WRITEdiagfi(ngrid,"raveco2","ice eff radius","m",0
2876	& ,raveco2)
2877
2878	endif ! of if (co2clouds)
2879
2880	CALL WRITEDIAGFI(ngrid,'reffice',
2881	& 'reffice',
2882	& 'm',0,rave)
2883
2884	do iq=1,nq
2885	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
2886	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
2887	end do
2888
2889	call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
2890	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
2891	call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
2892	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
2893	call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
2894	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
2895	& +zdqcloud(1,:,igcm_h2o_vap))
2896
2897	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
2898	& rice)
2899	ENDIF ! of IF (water)
2900
2901	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2902	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2903
2904
2905	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
2906
2907	do l=2,nlayer-1
2908	tmean=zt(1,l)
2909	if(zt(1,l).ne.zt(1,l-1))
2910	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
2911	zlocal(l)= zlocal(l-1)
2912	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
2913	enddo
2914	zlocal(nlayer)= zlocal(nlayer-1)-
2915	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
2916	& rnew(1,nlayer)*tmean/g
2917	#endif
2918
2919	END IF ! if(ngrid.ne.1)
2920
2921	icount=icount+1
2922
2923	END SUBROUTINE physiq
2924
2925	END MODULE physiq_mod

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

Download in other formats:

Original Format