Context Navigation

physiq.F @ 1114

Last change on this file since 1114 was 1114, checked in by emillour, 11 years ago

Mars GCM:

Correction: in physiq.F, move newcondens (CO2 condensation routine) so that it is the last atmospheric physical process to be computed.
Also corrected the setting of the South Pole surface temperature to CO2 condensation temperature (when obliquity < 27) to account for varying CO2 mixing ratio. This operation is now done in newcondens.

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

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

Download in other formats:

Original Format