Context Navigation

physiq.F @ 833

Last change on this file since 833 was 833, checked in by jbmadeleine, 12 years ago

Mars GCM:

Implemented the thermal inertia feedback:

Added a flag in callphys.def called tifeedback, set to false by default;
Changed physiq.F to call soil.F with a new thermal inertia if tifeedback = true;
Added a routine called soil_tifeedback.F that computes the new thermal inertia of the subsurface when ice is deposited;

Modified files: soil.F, physiq.F, inifis.F, callkeys.h
Added files: soil_tifeedback.F

File size: 91.0 KB

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

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

Download in other formats:

Original Format