Context Navigation

physiq.F @ 1502

Last change on this file since 1502 was 1502, checked in by emillour, 9 years ago

Mars GCM:

added missing condition to use MY32 dust scenario in aeropacity
corrected stats output for vmr_ice*rice in physiq

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

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

Download in other formats:

Original Format