Context Navigation

physiq.F @ 705

Last change on this file since 705 was 705, checked in by emillour, 13 years ago

Mars GCM:

Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

FGG

File size: 81.9 KB

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

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

Download in other formats:

Original Format