Context Navigation

physiq.F @ 722

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

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

Download in other formats:

Original Format