Context Navigation

physiq.F @ 633

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

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

Download in other formats:

Original Format