Context Navigation

physiq.F @ 370

Last change on this file since 370 was 358, checked in by aslmd, 14 years ago

LMDZ.MARS:

A FIRST VERSION WITH SCAVENGING
TRANSPARENT TO CASUAL USER
FOR WATER CYCLE FOLKS, SEE BELOW

[TN and AS solved conflicts due to r330 >> r357]

07/11/11 == JBM

Changed watercloud.F to call two separate routines,

simpleclouds.F or improvedclouds.F, which are a simplified and
full microphysical scheme for cloud formation, respectively.
Removed the tag called "improved" in watercloud.F, and added
another tag called "microphys" which is defined in callphys.F
instead. Changed routines: callkeys, inifis, physiq, watercloud.

Reimplemented the use of typical profiles of dust particle sizes

and CCNs in simpleclouds.F. Corrected the previously used
analytical CCN profile. Moved ccn_factor to simpleclouds.F,
which won't be used in the new microphysical scheme. Changed
routines: aeropacity, physiq, simpleclouds, watercloud.

Computed rdust at the end of callsedim instead of updatereffrad,

except at firstcall. Renamed rfall into rsedcloud and moved it
in simpleclouds. Moved nuice_sed to initracer.F and added it to
tracer.h. Changed routines: callsedim, physiq, tracer.h,
watercloud, initracer, simpleclouds, updatereffrad.

Added two tracers for the CCNs. Added the different tests in

initracer.F to make sure that, for example, microphys is not
used without doubleq, etc. Corrected an inconsistency in
callsedim.F, and changed the way r0 is updated. Changes
routines: callsedim, inifis, initracer, physiq, testphys1d,
tracer.h.

Added the ability to have a spatially variable density in

newsedim (same method as that used for the radius of
sedimentation). Required the addition of one input to newsedim,
which is the size of the density variable "rho". Changed
routines: callsedim, newsedim.

Added an output to aeropacity called "tauscaling", which is a

factor to convert dust_mass and dust_number to true quantities,
based on tauvis. Removed ccn and qdust from aeropacity, which
are obsolete now.

Wrote improvedcloud.F which includes all the microphysical

scheme, and connected it to the sedimentation scheme. Added and
changed routines: callsedim, physiq, growthrate, nuclea,
improvedclouds, initracer, watercloud, microphys.h.

07/11/11 == TN

Added CCN & dust tracers updates in physiq.F

Corrected a bug that can give negative CCN numbers, due to the
use of single precision values (only 6 decimals) whereas up to 10E+10
CCN can be formed or disappears...

Corrected physical bug that causes h2o_vap or h2o_ice

to be negative in improvedclouds.F.

Corrected physical bug that causes CCN & dust tracers

to be negative in improvedclouds.F when all ice sublimates and
releases dust

Added parameter contact mteta in callphys.def

Default value is still 0.95, see inifis.F

Changed tendancy computation for dust_number in improvedclouds.F

that was not the right one. Indeed, the scavenged dust_number tracer
is computed from the dust_mass one, and its tendancy before scavenging
must be taken into account to compute its scavenging's tendancy.

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

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

Download in other formats:

Original Format