Context Navigation

physiq.F @ 672

Last change on this file since 672 was 669, checked in by tnavarro, 13 years ago

Ice effective radius for output is weighted by ice total surface instead of ice total mass in physiq.F

New distribution for permanent ice reservoirs in surfini.F

icelocationmode = 1 allows for realistic ice distribution, no matter what resolution.
icelocationmode = 2 is predefined 32x24 or 64x48 and should be USED BY DEFAULT. It currently overestimates ice.
icelocationmode = 3 is good for quick logical definitions (paleoclimates,stability studies,etc...)

Possibility to change perihelion date in solar longitude as well as in sol in tabfi.F for newstarts
Improved latent heat release when ground ice sublimates: now a implicit scheme for ground temperature

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

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

Download in other formats:

Original Format