Context Navigation

physiq.F @ 682

Last change on this file since 682 was 674, checked in by acolaitis, 13 years ago

LES RESTART
===========

Several corrections for LES restart. Added a number of save variables from physiq to restart netcdf, so that continuity between runs can truly be achieved.
Added some lines in makemeso so that debug option "-g" also works for WRF (in addtion to the GCM physics) when compiling with ifort

===============================
WARNING WARNING WARNING WARNING
===============================

FROM THIS REVISION, YOU MUST MODIFY MANUALLY SOME FILES AS FOLLOW, BEFORE YOU CAN RE COMPILE THE LES :

simply copy the call_meso_physiq*.inc files in $MMM/SRC/WRFV2/ into $MMM/SRC/LES/WRFV2/.
=>> this is usually done during LES installation by $MMM/SRC/LES/LMD_LES_MARS_install

IF YOU WANT TO MAKE A RUN WITH RESTART FILES, YOU MUST MAKE A CLEAN RECOMPILATION OF THE MODEL WITH THE UPDATED REGISTRY AND UPDATED GCM FILES

=>> it is advised to remove your working directory (ex: lesnewphys_mpifort64) and start again with makemeso

IF YOU WANT TO MAKE A RUN AND GENERATE RESTART FILES, YOU MUST RECOMPILE THE IDEAL.EXE STEP AND PERFORM IT

=>> current wrfinput will not work, as it does not contain the new variables
===============================
WARNING WARNING WARNING WARNING
===============================

File size: 80.7 KB

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

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

Download in other formats:

Original Format