Context Navigation

physiq.F @ 1268

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

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

Download in other formats:

Original Format