Context Navigation

physiq.F @ 1226

Last change on this file since 1226 was 1226, checked in by aslmd, 11 years ago
LMDZ.MARS : Replaced comcstfi and planete includes by modules.
File size: 92.1 KB

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

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

Download in other formats:

Original Format