Context Navigation

physiq.F @ 1233

Last change on this file since 1233 was 1233, checked in by aslmd, 11 years ago
LMDZ.MARS. Filling geom arrays is now out of phys_var_state_init. Done through a merged function ini_fillgeom within the comgeomfi_h module. Cosmetic changes. New interface with the mesoscale model: lesser amount of dirty MESOSCALE includes.
File size: 92.9 KB

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

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

Download in other formats:

Original Format