Context Navigation

physiq.F @ 1528

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

Mars GCM:

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (==jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Added "ioipsl_getin_p_mod.F90" (getin_p routine) in phy_common to correctly read in parameters from *.def files in a parallel environment.

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

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

Download in other formats:

Original Format