Context Navigation

physiq.F @ 1467

Last change on this file since 1467 was 1467, checked in by tnavarro, 9 years ago
new option supersat to allow for supersaturation of water
File size: 100.3 KB

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

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

Download in other formats:

Original Format