Context Navigation

physiq.F @ 1377

Last change on this file since 1377 was 1377, checked in by emillour, 10 years ago

Mars GCM:

Update pbl_parameter with corrected version by Anne-Flore Moreau.

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

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

Download in other formats:

Original Format